{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <div align=\"center\">一元线性回归(SLR)</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "\n",
    "from scipy import stats\n",
    "\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings('ignore')\n",
    "    \n",
    "%matplotlib inline\n",
    "\n",
    "plt.style.use('ggplot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. 读取数据\n",
    "\n",
    "贾俊平<<统计学(第六版)>>线性回归章节例题数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>不良贷款\n",
       "(亿元）</th>\n",
       "      <th>各项贷款余额\n",
       "(亿元)</th>\n",
       "      <th>本年累计应收贷款\n",
       "(亿元)</th>\n",
       "      <th>贷款项目个数\n",
       "(个)</th>\n",
       "      <th>本年固定资产投资额\n",
       "(亿元)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.9</td>\n",
       "      <td>67.3</td>\n",
       "      <td>6.8</td>\n",
       "      <td>5</td>\n",
       "      <td>51.9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.1</td>\n",
       "      <td>111.3</td>\n",
       "      <td>19.8</td>\n",
       "      <td>16</td>\n",
       "      <td>90.9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.8</td>\n",
       "      <td>173.0</td>\n",
       "      <td>7.7</td>\n",
       "      <td>17</td>\n",
       "      <td>73.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3.2</td>\n",
       "      <td>80.8</td>\n",
       "      <td>7.2</td>\n",
       "      <td>10</td>\n",
       "      <td>14.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>7.8</td>\n",
       "      <td>199.7</td>\n",
       "      <td>16.5</td>\n",
       "      <td>19</td>\n",
       "      <td>63.2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   不良贷款\\n(亿元）  各项贷款余额\\n(亿元)  本年累计应收贷款\\n(亿元)  贷款项目个数\\n(个)  本年固定资产投资额\\n(亿元)\n",
       "0         0.9          67.3             6.8            5             51.9\n",
       "1         1.1         111.3            19.8           16             90.9\n",
       "2         4.8         173.0             7.7           17             73.7\n",
       "3         3.2          80.8             7.2           10             14.5\n",
       "4         7.8         199.7            16.5           19             63.2"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "db = '/home/lidong/Datasets/'\n",
    "bad_loans_df = pd.read_excel(os.path.join(db, \"Statistics/bad-loans.xls\"), usecols=\"B:F\")\n",
    "bad_loans_df[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 不良贷款Y, 各项贷款余额X\n",
    "ys_loans_bad = bad_loans_df.iloc[:, 0]\n",
    "xs_loans_balance = bad_loans_df.iloc[:, 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. 画出散点图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f4d5bfaf438>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAEkCAYAAADU9rPiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcZFV5//HPqe6enmZfGojNiDqiJDL+QpRRYtxAFAVEFPIEEySoMBLUgHEnqIMiIIgRXKJj+DlRFPOgQhaC8FMUF0QWUWACShwMwdFAD6CQKXq6p87vj3Nrpqanl+ql+ta99/t+vfrVVffeqnueujP91Dn3LCHGiIiIiBRTLe8CiIiIyOwpkYuIiBSYErmIiEiBKZGLiIgUmBK5iIhIgSmRi4iIFJgSuUiHhRD+TwghtHHcHiGEWvb4WdnvoXZem4cQwl7Z794QwjOzx3vnWyqR6lEiF+m8DwPHTHVACGFP4Dpg/2zT5SGERcDrgKtCCAPjjl8WQriy5fnOIYQb57fYU5b36cD1IYQnkP6ONMvytyGEz4//8hFCODKE8PGW588MIXgb59kphPCzEMK10xx3Sghh2cwjESm+oAlhRCYXQjgLePMEu/qBUaDRsu2yGONbQwhvAd7Ssn0XoA6MtGw7IcZ4U3aOQeBq4O9ijF/Otn0F+Ffg58BLgf8GhklfCA4AFgN7A7/I3u9K4LnA21rOsT7G+PBMY55OCGEp8P+yGH6QJe0bgDOyOF8K3Af8Z7ZtL2BnYKcsDoB/B/YEPtry1r+JMT7Wcp7FwBWkz+D3gRtjjB8YV5aVwKHAYHaOXwBrSJ8FwNOBtcAY8PUY4wfn/gmIdJfevAsg0s2yxDE+ebwccODpMcbhCV42CHwyxvjJid4zhHApsF32+DnAF0nJ/sEQwo+BCGwALgE+T0rgTwD6Y4wnZa9bBpwdYzw6e34usANwdnaaZwCXAufPLvKJhRAOBy4GdgUWt5R3E/BV4DJgPfBE4OEY4xHZ644EDo0xnp49vwzYsaW8y4H3kj5XQgi7ApcDD5O+nGwPfDuEsCPwzhjjpux1TwY+BmwXY7w0hLAzKfm/CFgE/Ag4OCvjI/P5WYh0CyVykRkIIewAfBb4yCRJfKYOAAy4ELgrxvislnN9BthISnb7AeeGEP4HuIDUnL0ohHA3qUn+EFIy/GqM8bAQwjdJtfz5tj/wEuCaGOO3gOa9/ABcRaqJXwgcARwWQhgCTgd6gJ7sS9CVWdwHAxfEGF8XQvgP4FvZe72A9OXm34G3xhgbwKMhhINJyf2GEMKKGONPY4wnhhC+AewbQjglK+Mbs7JsDzyUfQ77AK8Bvt+Bz0QkV0rkIm3KktUXSE3F10xz+BktiWW8JcA/AMQYV2Xv/VO2bnoHOBVYCTT3fYLUpP2WccdtBO4HHiMltB7g1hjjHdNHNTMxxguy8t40bnsMIbyGVDv+CenWw6nAraQm9la9wL2k2wNLstr3l2KM60MIFwOvICX1w4HDQwi7kG5h/C57/eeBfwshfJj0BeFh0peLPYFPZu8N6YvCy0lfdN4/Lx+ASBdSIhdp37nAc4CvAUeGEE4GfgZc1NLU23TONE3r4y0jNR1P9JJFwPtjjB5C2B74vXH7fxZjPCdLePsAi2KM754umBDCq4BzgG/EGN+eJdS3xRjbSXrPCiHcOcm+nYFjYow3ZT3bx5f31hjjN0IIB5BuTzxM6hAIqdn+b2KMYy3lPBsYjjG2dpb7CKlV4sjstf9CamZ/Q4xxY8vneCzpi5BIaSmRi7QhhHA+8Jeke6/HZ5vPBVYBx4UQToox3t7ykg+EEE4n3QvfEfifln13An3jTrFfjPHJ2bn2AZ6WNV0TQlgNPJAddz5bOnK1lu9EUgIdJTWxf7eNsM4jdRT7YAjhvaRm53Z7fu/QUt5lwECM8ebs+XdayvsZ0v3y8eU9jNRMv2sIYZ8Y433ZrlcDp437QrMT0AghvGPc2xwHPAqcBdwCvBX4mxDCi1uO6QcebzMmkULS8DORKYQQdsx6kL8aeEmM8e7mvhjjL2OMLwNWA98LIbTWgs+KMe5L6jV9L2DZ89eTkvhfZr2yJ7IMOHqiHTHGN5Oa2F9P6gT3/qxcd5KapM8B/qLN8K6NMf6K1AT+h6Rm6dn06n5h9jNReV9NSuZnkpLtR4CDYozrJypvjPGCGOOSrDx/lD2+GFiZPX5xjHFJ9vN90m2KPyB9KTifVPs/o6UIe5I634mUlhK5yNQ+Q2q5enaWLLcRY/w08Dzg9gn2PUYavnZ11sP7LcC7Y4zHxxgnqykuJ/UI7x+/I2uOPp5U8wbYA/hQCGF34DBSojwgq9VPKcZ4WvZ7JMZ4XIzxD2OMP5nudZOUd4fs3vz48h4FHNSy6RnAX4cQnkbqif8xUovGTuNeugrY6vZACGFf4IchhENaYvgc8FCM8aAY40Gk2vcoqUVgZ2Bf0henh9m2D4JIKSiRi0ztr2KMx8YYfzfVQTHGNTHGZi/xPmDHEMIrQwifInW4ejmpif2pwLLsfnar40IIu4cQzgBeSxp+dmcI4bOkYV0/yo47hdQLfEfSkK8vAt8hdcI7N8Y4ShqbvXqKGv98OC6EsHcI4UJSbXgf4K7s+dmkzneQbkf8bUt5LySN674UeG+McUP2+NNhy6x2pwJPA97XesIY43+SJsi5Imuan8qfkT6HHwDXk/oB3Dy3kEW6VIxRP/rRzwx+SIlq5RT7byHVAC8D/hzYNdveC7yBlFxGSMPI+kjjxX9M6jh3FrBLdnw/KSFdn+3fI9v+NVKN89Ts+eeB88eV4XPAxzsQ+67AV0j3+W8H/oZ0fxzSOPaTSD3VryN1ugO4MSvvUdnz64A3t7xnD/AN0heUo4FfA09t2f9Z4KSW50eRerA/M3v+y+wcN2bbD8kefwgIpOb3a4A35f1vRz/66cSPZnYTmaGsF/VYjHHlJPufCdwTJ286J5vYpDfG+HAI4UDggbilw9dExz8lxnjvJPueEGP89bhtNdIkKY9N9Jq5CCH8CfCLGONvpjhmaYxx7ST7JirvYtIQsyXAbjHGW7IpXF9AGgv+ktbzhRD+GPhRjLERQtgtxvhQy75nABtjqsE3t9WAvWOMzZnlREpDiVxERKTAdI9cRESkwJTIRURECkyJXEREpMCUyEVERApMiVxERKTAlMhFREQKTIlcRESkwJTIRURECkyJXEREpMCUyEVERApMiVxERKTAlMhFREQKTIlcRESkwJTIRURECqw37wLMgNZbFRGRqgnTHVCkRM66desAGBwcZHh4OOfSdFYVYoRqxFmFGKEacSrG8ihCnENDQ20dp6Z1ERGRAlMiFxERKTAlchERkQJTIhcRESkwJXIREZECK1SvdRGZX416g01rR4l1CAPQs7SP2oC+34sUiRK5SEU16g3Gbh4hBgi1QKMeietH6F3er2QuUiD63ypSUZvWjm5O4pB+x5C2i0hxKJGLVFSsb0niTaEWiPWcCiQis6JELlJRYQBiY+uZj2MjEgZyKpCIzIoSuUhF9SztI8QtyTw2IiGm7SJSHOrsJlJRtYEavcv7N/darw0E9VoXKSAlcpEKqw3UqO3fn3cxRGQO9NVbRESkwJTIRURECqzjTetmthg4GXg7sNLdV5vZkcA7gB2Bq939zE6XQ0REpIwW4h75XkAd+FLLttuBlwIN4F4zu9DdHx7/QjNbAawAcHcGBwcB6O3t3fy4rKoQI1QjzirECNWIUzGWR5niDDHG6Y+aB2a2Evilu69u2bYLcBOwv7tPN51UXLduHQCDg4MMDw93qKTdoQoxQjXirEKMUI04FWN5FCHOoaEhgDDdcbndIzez7YEvA6e2kcRFRERkArkkcjMbBK4APuLu38yjDCIiImWwEJ3dlgBXAkPAiJkdAuwKPA04y8wAVrn7lztdFhERkfmW93LAHU/k7n4/cGCnzyMiIrLQumE5YI0jFxERmaVuWA5YiVxERGSWumE5YCVyERGRWeqG5YCVyEVERGapG5YD1upnIiIis9QNywErkYuIiMxB3ssBq2ldRESkwJTIRURECkyJXEREpMCUyEVERApMiVxERKTAlMhFREQKTIlcRESkwJTIRURECkyJXEREpMCUyEVERApMiVxERKTAlMhFREQKTIlcRESkwJTIRURECkzLmIqIdFCj3ti8VnUYYMHXqpbyUyIXEemQRr3B2M0jxAChFmjUI3H9CL3L+5XMZd7oX5KISIdsWju6OYlD+h1D2i4yX5TIRUQ6JNa3JPGmUAvEek4FklJSIhcR6ZAwALERt9oWG5EwkFOBpJQ6fo/czBYDJwNvB1a6+2ozew7wUaAH+KG7v6PT5RARWWg9S/uI60dS8q6F9Dum7SLzZSE6u+0F1IEvtWxbBRzl7veZ2TfN7HnufsP4F5rZCmAFgLszODgIQG9v7+bHZVWFGKEacVYhRqhGnLOJcdPuY9Tv3kDjfxvUtq8x8Pvb0bNd9/YzrsJ1hHLFGWKM0x81D8xsJfBL4J+Ba4HDgAuA/YHL3P2iad4irlu3DoDBwUGGh4c7VtZuUIUYoRpxViFGqEacirE8ihDn0NAQQJjuuLy+Fg4BlwLvBV6I7tWLiIjMyoInUHd/GHgMOAm4HTgc2KZZXURERKa3EJ3dlgBXkmrhI2Z2CPBXwBXAKHCtu9/c6XKIiIiUUccTubvfDxw4wa7ndvrcIiIiZad70yIiIgWmRC4iIlJgSuQiIiIFpkQuIiJSYErkIiIiBaZELiIiUmBK5CIiIgWmRC4iIlJgSuQiIiIFpkQuIiJSYErkIiIiBaZELiIiUmBK5CIiIgWmRC4iIlJgSuQiIiIFpkQuIiJSYErkIiIiBaZELiIiUmBK5CIiIgWmRC4iIlJgSuQiIiIFpkQuIiJSYErkIiIiBaZELiIiUmBK5CIiIgXWm9eJzWwx8AVgCFgMnOPuX8+rPCIiIkWUWyIHnguMuPvzzezZwAeArRK5ma0AVgC4O4ODgwD09vZuflxWVYgRqhFnFWKEasSpGMujTHGGGGMuJzazfuArwHrgScDfuvtNU7wkrlu3DoDBwUGGh4c7X8gcVSFGqEacZY6xUW+wae0osQ477bkjG/Z4nNpAee/YlflaNlUhRihGnENDQwBhuuPy/B+3E7A9cBMwCizPsSwiMkONeoOxm0doPNSAkcjog6Ppeb2Rd9FEKiXPRH4UcLu7ryI1n5+cY1lEZIY2rR0lBgi1VGEItUAMabuILJw875H/O/AXZvadrBwfyrEsIjJDsb4liTeFWiDWcyqQSEXllsjd/dfAIXmdX0TmJgxAox63SuaxEakNTHtLT0Tm0bSJ3Mx2AU4Crnb3NRPsXwScB5zp7hvmv4gi0o16lvYR148QGymZx0YkxLRdRBZOO/fIdwdeB3zczM5t3WFmS4Crge2VxEWqpTZQo3d5P7XdatAf6NujLz0vca91kW7UbtP6f7j7a80sZMPGlgPHA38IXOjuX+1YCUWka9UGatT27wdgh8GdeHx4Y84lEqmetr86m9nVwD3Ao8A3gPuA5yuJi4iI5GfKRG5mPwQ+BuwJHOvu+wJ7AC8DBoAbzOxVHS+liHSFRr3B6JoRNt4ywugajRkX6QbT1cgPI03Y8iPgZjN7j7v/FnjE3d8HvAI42czO63A5RSRn4yeAaTzU0AQwIl1gykTu7r8DXuXuZwDfBr5vZgcCq81sR+D9wCXAnma2V8dLKyK50QQwIt1pJuPIDwa+RpoXPQB/Atzq7lcAV3SgbCLSRTQBjEh3mu4e+RXAU83s66SlRludBvy5mf2rmb2pUwUUke4QBtKEL61iIxIGciqQiADT3yN/H/Cr7Hdr+1kE3kFK5h8EzupI6USka/Qs7SPELclcE8CIdIfp7pHfCTyezeg2vkfLs4FPAI8B3+9M8USkW4yfAKa2W00TwIh0gXbukf80+/1OUsLuBf7Y3b9gZrcDb3X3YztVQJFu0br2dhhINdGqJbHWCWBEpDu0k8i/ZmZvAF4N7OHunwfeAuDuPwFONbPd3X19B8spkqvm0Ktmr+1GPRLXj6hGKiK5mzKRZ9OxDgK7AWcAXzCzR4EHWg67CbjGzI5x9//qWElFcjTh0KtGZNPaUdVQRSRX09XIvwb8gNRj/e7s93nA97L9EdgBeExJXMpMQ69EpFtNmsjN7AnAd4HrgNcD5wMOvNrdX99y3N8Dn+pwOUVypbW3RaRbTXpzz91/TZpn/enAG4B+YCVpMhhg81rkLwGu7GgpRXKmoVci0q2mqpH/EXA5cDNwk7ufmm3HzA4mdX7rB65yd83RKKXWHHrV7LVeGwiV7LUuIt1n0kTu7rcB+wKY2c3jdt8DXE1al/w4M9spm5ddpLQ09EpEutFUNfIAPCObDKa12hHd/X7gfuBqMxsF3oZmdxMREVlwU/Va3x24zMzOAd5lZkPAxgmO+yypZ7sSuYiIyAKbqml92MwOBf4NeCJwB3A7sM7M7s0Oi8BrgQ1mtiSrqYuIiMgCmXIcubs/YGaHk6Zmfbu73zHRcWZ2G7AnqbldREREFsi0XW7dfRh4HfDcKQ47xd1/PG+lEhERkba0NXbG3W92939o3ZaNIW/6sJntPK8lExERkWlNu2iKmR0IbDducx/wl8AJWQL/a9K65DNiZs8hdZLrB77k7pfM9D1ERESqLMQYpzzAzO4iTQoTgCNInd92ICXzdwPPJC1l+vyZnNjMdgSuBY5y9wcnOWYFsALA3Z+9cWPqNN/b28vY2NhMTlc4VYgRqhFnFWKEasSpGMujCHEuWrQIWmZTnUw7y5ji7icAmNnt7n6Cmd1Dmn/9/cBS4NOzKOMLgEXAP2ZJ/QJ3/5dx510FrMqexuHhYQAGBwdpPi6rKsQI1YizCjFCNeJUjOVRhDiHhobaOq6de+QTVdmju38feBqwM/CF9ou22U7ALe5+OPAa4KJZvIeIiEilTZnIzewY0v3rifYdAOwNmLs3ZnHue4AnZY83Av87i/cQERGptEkTuZk9k9ShbWCC3b8HrAaOBY6dTY91d78VuNvMbgD+BThtpu8hIiJSdVPN7HYHcFTW2W28R4HnufsGMzsOeAXwlZme3N1Pn+lrREREZIvZ3iPf4O4bssffAF45f0USERGRdrXTaz2Y2amkLvC7Z493zn5D6rT2cjNb5O4TLaoiIiIiHdJOIv8ksEf2eFX2+FMt2/qAtaQe7Gvmu4AiIiIyuXYS+dPdfZuOaGb2TuB2d78GOHPeSyYiIiLTauce+WvHb8h6tJ8B3DfvJRIREZG2tZPIt5oezsxqwOeAM9x9oh7tIiIiskBm02v9NOC37v73HSiPiIiIzEA7ibzHzHZteb4fcHKHyiMiIiIz0E5ntz7gXjM7j7Ra2a+BE81smwPdfcZLmYqIiMjstVMjHwH+DDgF+D6wLym5T/QjIiIiC6idGnl092vMbBlplbMXAR9y9593tmgiIiIynXZq5AC4+2Pu/hrgy8D1ZvYHnSuWiMxEo95gdM0IG28ZYXTNCI36bBYkFJEiajuRN7n7e4AvAlea2Q7zXyQRmYlGvcHYzSM0HmrASKTxUPZcyVykEmY8jhzA3d8F3Ab833kvkYjMyKa1o8QAoZb+q4ZaIIa0XUTKr51E/qlJtp8EbFStXCRfsb4liTeFWiDWcyqQiCyoaTu7ufvKSbY/Bhw/3wUSkZkJA9Cox62SeWxEagPbNKaJSAnN+B65iHSXnqV9hJiSN6TfIabtIlJ+7Qw/E5EuVhuo0bu8P90rr0NtINCztI/agL6ni1SBErlICdQGatT278+7GCKSA31lFxERKTAlchERkQJTIhcRESkwJXIREZECU2c3kUyj3tjc8zsMoJ7fIlIISuQibJmvvDnVaaMeietH6F3er2QuIl0t10RuZm8BPg7s6+6/zLMsUm0TzlfeiGxaO6phXSLS1XJL5Gb2FOCVwA1THLMCWAHg7gwODgLQ29u7+XFZVSFG6J44f1d7hMYO264WVqvV2Glwlzm9d7fE2GlViFMxlkeZ4swlkZtZAD4BvBVYNdlx7r6qZX8cHh4GYHBwkObjsqpCjNA9cY42Rmg81th2vvJFNTYOj83pvbslxk6rQpyKsTyKEOfQ0FBbx+V1828FcJ27/zyn84tsRfOVi0hR5ZXIXwkcbWbfAQ4AvmJm++VUFpHN85XXdqtBf6C2W00d3USkEHJpWnf3I5uPs2R+ojq7Sd40X7mIFFHuw8/c/cV5l0GkDDQOXqSack/kIjJ3GgcvUl36Hy5SAhOOgw9pu4iUmxK5SAnEOlsNnYMsmddzKpCILBglcpESCANbhs41xUYkDORUIBFZMLpHLlPqdAeqTRvGGF0zog5ac9SztI+4fiQl72x6WY2DF6kGJXKZVKc7UDXqDX538yM0Hm+og9YcNcfBN7901QaCvhSJVIQSuUyq0wuJbFo7Sl9tsRYqmScaBy9STUrkMqFGvcHYvWNQj9ALYdcatb7avHaginUIPeqgJSIyF2p3k21sblLf2CCORmI9En+1icZoY147UKmDlojI3CmRyzaaTeo9u/XSrC/HAPGhxrx2oOpZ2geNqIVKRETmQIlcttEckxz6ArWhXsJAesyiMK8d0WoDNXZ60S5aqEREZA50j1y2EQagUY+bk3nPnr1pbe7davOeZHu266VPHbRERGatkolci0tMTWOSRUSKo3KJXItLTE9jkkVEiqNyibzTY6OLYrpWCY1JFhEphsolci0uoVYJEZEyqdxfbY1d1pKXIiJlUrlE3rO0jxCp9NhltUqIiJRH5RJ5syNXlccuq1VCRKQ8KnePHNSRS8PLRETKo5KJvOo0vExEpDyUyCuq6q0SIiJloSqYiIhIgSmRi4iIFJgSuYiISIHldo/czLYDVgF7A9sBJ7j7z/Iqj4iISBHlViN39w3AB9z9YOAS4NS8yiIiIlJUufZad/dfZA+HgLXj95vZCmBFdiyDg4MA9Pb2bn5cVlWIEaoRZxVihGrEqRjLo0xxhhjj9Ed1kJkdDbweOMbdx6Y4NK5btw6AwcFBhoeHF6J4ualCjFCNOKsQI1QjTsVYHkWIc2hoCCBMd1yuNXIzeyNwCPCn0yRxERERmUBu98jN7EBSZ7cnAtea2bV5lUU6q1FvMLpmhI23jDC6ZoRGvdGR14iIVFFuNXJ3vwXoyev8sjBms/a51ksXEWmf/ipKR81m7XOtly4i0j4lcumo2ax9rvXSRUTap0QuHTWbtc+1XrqISPuUyKWjepb2EeKWxNzO2uezeY2ISFVpGVPpqNmsfT7+NfREoMbYmlHCAFo7XUSkhRK5dNxs1j5vvmZLD/aGerCLiExAfwmlq7X2YI+jkcbwJsZ+M8bG79U1tlxEBCVy6XLNHuxxNNJYN0asR8KmQHwk1dSVzEWk6tS0PoVGvbH5Pq3uzc7eXD7HMACNeqTx8CYiEEIgxkjo2zK2fKbN9iIiZaJEPomyzC6W95eRuX6OPUv7iOtHiKNxSxKPEHataWy5iAhqWp9UGWYXaybRxkMNGIk0Hlr45ui5fo7NHuy1XWrQEwkDgbB3D7W+msaWi4igGvmkun12sXZq2hMm0UZc0Obo+fgcawM1Fr1gYKuavcaWi4gkqpFPoptnF2u3pt0NX0bm63PcXDPfrQb9gdputcLd5hAR6QTVyCex+d5sI3ZdDbDdmnazo1hrMo+NSG1g2nXq5818fo6zGY8uIlJ2qs5MoptrgO3WtLthqtNu/hxFRMpANfIpdGsNsN2a9mymR+2Ebv0cRUTKQIm8gGbSXK0kKiJSbkrkBTRfNe28x5iLiMjcKZEX1Fxr2mWZ8EZEpOqUyEuoKGPMRURk7iqXyMvenNxuTbsbxpiLiMjclSeDtaEbpizttHanRO3mCW9ERKR9lUrkZZg/fTpFGmMuIiJzV6lEXoXm5HZr2pqoRUSkHCp1j7wbpiztNI0xFxGpllwSuZn1ABcDy4Ae4K/c/Y5On7eb50+fL90ym5uIiCyMvP66HwvU3P1FwHuACxfipFVpTq4N1Ojbv59FB/bTt3/54hMRkS1CjHH6o+aZmV0EXEtWGwee4e5PmuC4FcAKAHd/9saNGwHo7e1lbGxs4QqcgyrECNWIswoxQjXiVIzlUYQ4Fy1aBDDtvd8875GfCVwNHAXcPdEB7r4KWJU9jcPDwwAMDg7SfFxWVYgRqhFnFWKEasSpGMujCHEODQ21dVxeba63Are5+weBA4G7ciqHiIhIoeVVI78MeJmZXQ+MAafkVA4REZFCyyWRu/socHwe5xYRESkTdWcWEREpsFx6rc9SYQoqIiIyT6bttV6kGnlo/pjZra3Py/hThRirEmcVYqxKnIqxPD8FinNaRUrkIiIiMo4SuYiISIEVNZGvmv6QwqtCjFCNOKsQI1QjTsVYHqWJs0id3URERGScotbIRUREBCVyERGRQlMiFxERKbA8Vz+bETPrAS4GlpEtf+rud+RbqvllZo8DN2ZPLwJ+BXyUFO8P3f0deZVttsxsMXAy8HZgpbuvNrPnMC6uIl/fSWJcCRwH/Ab4rbu/qsgxApjZdqQOQnsD2wEnADtTrms5UYyvpWTXEsDMrgEWAYtJy0UPUKJrCRPGeAwlvJZFqpEfC9Tc/UXAe4ALcy5PJ/zG3V+c/VxB+oNyvLv/CXCAmT0v5/LNxl5AHfhSy7aJ4iry9Z0oRoDzsmv5qux5kWPE3TcAH3D3g4FLgFMp2bWcJEYo2bUEcPfDsjivBp5Hya4lTBgjlPBaFqbXupldBFxL9o0JeIa7PynfUs0vM/tPYB0wDHwY+AxwGHABsD9wmbtflF8JZy+rof4S+GfSddwqLmApBb++zRizGvnppIWB6sDl7n5xmf4NZ7FG4EhKeC1hc4wPk2bXKt21NLMjgLOAXYAXAVdSsms5LsYXAkYJr2WRauQAZwIHAEeRlj8tFXff191fCFwKnA0MZY8vZtvaXpFNFldprq+7f9zdDwReChxvZvtluwofo5kdDfwR8GlKei1bYvxUWa+lu1+VxfXXpEpD6a7l+BjLei2LlMhvBW5z9w8CBwJ35VyeTgrAvcBjwEnA7cDhwA15Fmo+uPvDTBxXWa/63J0AAAAFWUlEQVTvpuz3o5QgRjN7I/CnwJ+6+4OU8FqOi7H1j3qprmWLR4ANlPBatniEVAtvKtW1LFLTeh/weeCJpG9Mp7j7PfmWav6Y2Z7A14AR4AHgzaQawbnAKHBt9g+tUMxsCanJbogU2/eA1YyLq8jXd5IY7wJeRmp+/oy7e5FjBDCzA4EfAT8AGsBG4DzKdS0nivHblO9aPoX0//Bx4HfA6cB+lOtaThTjCZTsWkKBErmIiIhsq0hN6yIiIjKOErmIiEiBKZGLiIgUmBK5iIhIgSmRi4iIFJgSuUiJmNkeZlbLHj8r+z1kZiHfkolIpyiRi5RENhfBdaTpNQEuN7NFwOuAq8xsoOXYZWZ2Zcvznc3sRqZhZn1m9j0zW5ONv53suKPN7NBZByMibdM4cpESMLNB0sIQf+fuX862fQX4V+DnpCkp/xtYSZo/fDFpha9fZG9xJfBc4G0tb7s+m4mveY4a8DnSvNWPkCoCb3T3Rssxp5Dmst6RNMHGfwDrgT7ShDlPJk14tAG40d1PmaePQKSyVCMXKbhsWdgfAPsAD5rZj83sVlKivgR4A2mJyicA52ZzTR8HXO/uB2bPtwd2IM3xfzbwddLSrM1zDABfJM3+dTzwJmBP4MvZ0p9NTyDNwf6P7r6buz8fWEJarOJQ4L9I038eCpzWgY9DpHJUIxcpODNbQZpW9ELgRHe/v2XfZ0jTxr6LlIRvA+4hfYlfRJq+8jrgEGA58FV3P8zMvgm8zd3vMLNlpEU01gJ/7u717L37SDX0PwZWuPv1Lec8krTaHaQFKU4iTTncnO96L+C97n7pvH8gIhXTm3cBRGRu3H0VgJn9lJS0W51Kak5v7jsNuHvcMRuB+0mLZuxrZj3ArVkSfxfwVlIT/eHAXWa2A+lLwEPZ6y8DVpnZ10ktA08GXpG93z+TFuEA+Avg90hfCI6aa9wikiiRi5THMuDbZjbRvkXA+4GrgFeP2/czdz/HzHYhNc8vcvd3Z/v+Cfiku29oHmxmJwHL3P30lm1nkprmX0hanGJVVp6z3f2hljIdCnxrTlGKyFaUyEXKYz93fzKAme0DPM3dv5U9X03qZHY+qVPbVszsRGBn0spXy4HvZrueB1ww7svBdkCfmR077m3eCdwHvAd4kFTrPtnMrgD+Nzumn9ScLyLzRJ3dRMppGXD0+I3u/mbgE8DrgWFSLf0l7n4nqTn8HFITePP4y9x9CfBU4AXZ43cBl2SPD3b3JdnPZaQa/R+QOtp9lHRf/t1ssSepF7uIzBMlcpFyWg4sNrP+1o1mdgCp1/lotmkP4ENmtjtwGPAR4ICsRt/qHOCCce+1M/BdM9sq8QOj7n6Qux9Euh9fI91Pj6Qa/hjwKFtq6SIyB0rkIuVxnJntbmZnAK8ljdW+08w+C3yV1LP9FNI97B2BTaQhZd8BvkAamjZKqkmvNrPFAGZ2JPBnpI5zm7n7b0kd4C42szdMU7bTgBOBnwGXA//j7lfMNWARUSIXKbxstrVLSOO3byDdhz7I3U8jNbFfR7p//QPgfcAHSZ3ersomc3klsMbdvw7g7peTJoo5z8yeC6wGjnX3B7JT7kHqkY6730ZK5hea2cuy/T1mdmM2U9wRpJnmriF1hjsCeD7pS8eHOvSRiFSKOruJFJy7j5rZ3wMPuPt94/aNkHqe/5OZPcXdHwSOGfcWZ7j7r8dtexOpU9tiwNz9RjP7KOn++eOk5N88x4/M7CXAXdmmI9y9OTQNM3si8At3/0m2aYSUyJ88+6hFpEkTwoiIiBSYmtZFREQKTIlcRESkwJTIRURECkyJXEREpMCUyEVERArs/wOLqXB3IYrynAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,4))\n",
    "fig.suptitle(\"不良贷款 ~ 贷款余额\")\n",
    "ax.set_xlabel(\"贷款余额\")\n",
    "ax.set_ylabel(\"不良贷款\")\n",
    "ax.grid(True)\n",
    "ax.scatter(xs_loans_balance, ys_loans_bad, alpha=0.5, color='orchid')\n",
    "\n",
    "# fig.tight_layout(pad=1)  # Adjust subplot parameters to give specified padding."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. 回归分析 \n",
    "\n",
    "Model: y ~ x + c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>const</th>\n",
       "      <th>各项贷款余额\n",
       "(亿元)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.0</td>\n",
       "      <td>67.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>111.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>173.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.0</td>\n",
       "      <td>80.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.0</td>\n",
       "      <td>199.7</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   const  各项贷款余额\\n(亿元)\n",
       "0    1.0          67.3\n",
       "1    1.0         111.3\n",
       "2    1.0         173.0\n",
       "3    1.0          80.8\n",
       "4    1.0         199.7"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xs_loans_balance_with_intercept = sm.add_constant(xs_loans_balance)\n",
    "xs_loans_balance_with_intercept[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 拟合\n",
    "model_ols = sm.OLS(endog=ys_loans_bad, exog=xs_loans_balance_with_intercept)\n",
    "fitted_ols = model_ols.fit()\n",
    "\n",
    "# 预测\n",
    "xs_pred = np.linspace(start=xs_loans_balance.min(), stop=xs_loans_balance.max(), num=60)\n",
    "xs_pred_with_intercept = sm.add_constant(xs_pred)\n",
    "ys_pred = fitted_ols.predict(exog=xs_pred_with_intercept)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAEkCAYAAADU9rPiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4XFeB/vHvmRnNaDQjV9mJFVdZSiCNQBISAgRIJT1OnOOwCyw1S13YvrAsG3rb7C4suwtm+ZFlKfFJM5CQAgkJJYVUUkiR3GLHSWy5aorKzJzfH2eUKLJsy7ak0cy8n+fRoyl35p4z19Y7595TjPceERERqU6RShdARERE9p+CXEREpIopyEVERKqYglxERKSKKchFRESqmIJcRESkiinIRcaZMeZoY4wZxXazjDGR8u3XlH+3jua1lWCMOaj8O2aMOap8+5DKlkqk/ijIRcbfF4CL97SBMWY2cDtwRPmhq40xceAdwI3GmOSw7Y80xqwccn+qMeaesS32Hst7KHCnMWYO4e/IYFn+0RjzveFfPowx5xpj/n3I/aOMMW4U+5lijHnKGHPrXrb7gDHmyH2viUj1M5oQRmT3jDGfAT48wlMJYAAoDXnsx977jxpjPgJ8ZMjj04A80DfksXd6739f3kcLcBPwb977H5Ufuwr4GfA0cDqwHugmfCE4BmgEDgFWld9vJXAC8JdD9rHFe79tX+u8N8aYNuAX5Tr8rhzadwGfLNfzdOAZoKv82EHAVGBKuR4APwdmA/8y5K2f995nhuynEbie8Bm8ArjHe//Pw8pyOXAa0FLexyrgccJnAXAosBooANd57z974J+AyOQSq3QBRCazcnAMD4+3Ag441HvfPcLLWoBveu+/OdJ7GmN+ADSVb78W+D9C2G82xjwIeCAHfBf4HiHA5wAJ7/37yq87Evi89/7C8v0vAWng8+XdHA78APjq/tV8ZMaYs4FvANOBxiHlLQLXAD8GtgDzgG3e+3PKrzsXOM17//Hy/R8DzUPKezzwCcLnijFmOnA1sI3w5SQF/MoY0wz8rfe+WH7dQuBfgSbv/Q+MMVMJ4f8mIA7cC7ylXMbtY/lZiEwWCnKRfWCMSQPfBr6ymxDfV8cAFrgCeMJ7/5oh+/oW0E8Iu8OALxljXgC+RjidHTfGPEk4JX8KIQyv8d6faYz5JaGVP9aOAE4FbvHe3wYMXss3wI2ElvgVwDnAmcaYVuDjQBSIlr8ErSzX+y3A17z37zDG/BG4rfxebyR8ufk58FHvfQnoMca8hRDudxljLvPe/8F7/y5jzM1AuzHmA+UyvrdclhSwtfw5zAcuAn47Dp+JSEUpyEVGqRxW3yecKr5lL5t/ckiwDDcX+B8A7/3y8nv/gZefegf4EHA5MPjcfxBOaX9k2Hb9wAYgQwi0KPCA9/7Rvddq33jvv1Yu7++HPe6NMRcRWscPEy49fAh4gHCKfagYsIZweWBuufX9Q+/9FmPMN4CzCKF+NnC2MWYa4RLGzvLrvwfcYIz5AuELwjbCl4vZwDfL7w3hi8JbCV90Pj0mH4DIJKQgFxm9LwGvBa4FzjXGvB94Cvj6kFO9g764l1Prwx1JOHU80kviwKe9984YkwIOHvb8U977L5YDbz4Q997//d4qY4y5APgicLP3/q/LgfqX3vvRhN5rjDGP7ea5qcDF3vvfl3u2Dy/vA977m40xxxAuT2wjdAiEcNr+r7z3hSHl/DzQ7b0f2lnuK4SzEueWX/tTwmn293jv+4d8jksJX4REapaCXGQUjDFfBf6McO317eWHvwQsBy41xrzPe//IkJf8szHm44Rr4c3AC0OeewxoGLaLw7z3C8v7mg90lE9dY4y5EthU3u6rvNSRa2j53kUI0AHCKfZfj6JaXyZ0FPusMeYThNPOo+35nR5S3iOBpPf+vvL9O4aU91uE6+XDy3sm4TT9dGPMfO/9M+WnlgAfG/aFZgpQMsb8zbC3uRToAT4D3A98FPgrY8ybh2yTAHpHWSeRqqThZyJ7YIxpLvcgXwKc6r1/cvA57/1a7/0ZwJXAb4wxQ1vBn/HetxN6Ta8BbPn+uwkh/mflXtkjORK4cKQnvPcfJpxifzehE9yny+V6jHBK+ovAn46yerd6758lnAJ/FeG09P706j65/DNSeZcQwvxThLD9CnCi937LSOX13n/Nez+3XJ5Xl29/A7i8fPvN3vu55Z/fEi5TvJLwpeCrhNb/J4cUYTah851IzVKQi+zZtwhnro4th+UuvPf/BZwEPDLCcxnC8LWbyj28PwL8vff+7d773bUUjyf0CE8Mf6J8OvrthJY3wCzgc8aYmcCZhKA8ptyq3yPv/cfKv/u895d671/lvX94b6/bTXnT5Wvzw8t7PnDikIcOB/7CGNNB6In/r4QzGlOGvXQ58LLLA8aYduBuY8wpQ+rwHWCr9/5E7/2JhNb3AOGMwFSgnfDFaRu79kEQqQkKcpE9+6D3fqn3fueeNvLeP+69H+wl3gA0G2POM8b8J6HD1VsJp9gXA0eWr2cPdakxZqYx5pPA2wjDzx4zxnybMKzr3vJ2HyD0Am8mDPn6P+AOQie8L3nvBwhjs6/cQ4t/LFxqjDnEGHMFoTU8H3iifP/zhM53EC5H/OOQ8l5BGNf9A+AT3vtc+fZ/mZdmtfsQ0AH809Adeu+7CBPkXF8+Nb8nywifw++AOwn9AO47sCqLTFLee/3oRz/78EMIqsv38Pz9hBbgj4E/AaaXH48B7yGESx9hGFkDYbz4g4SOc58BppW3TxAC6c7y87PKj19LaHF+qHz/e8BXh5XhO8C/j0PdpwNXEa7zPwL8FeH6OIRx7O8j9FS/ndDpDuCecnnPL9+/HfjwkPeMAjcTvqBcCDwHLB7y/LeB9w25fz6hB/tR5ftry/u4p/z4KeXbnwMM4fT7LcCfV/rfjn70Mx4/mtlNZB+Ve1EXvPeX7+b5o4BOv/tT55QnNol577cZY44DNvmXOnyNtP0i7/2a3Tw3x3v/3LDHIoRJUjIjveZAGGNeD6zy3j+/h23avPerd/PcSOVtJAwxmwvM8N7fX57C9Y2EseCnDt2fMeZ1wL3e+5IxZob3fuuQ5w4H+n1owQ8+FgEO8d4PziwnUjMU5CIiIlVM18hFRESqmIJcRESkiinIRUREqpiCXEREpIopyEVERKqYglxERKSKKchFRESqmIJcRESkiinIRUREqpiCXEREpIopyEVERKqYglxERKSKKchFRESqmIJcRESkisUqXYB9oPVWRUSk3pi9bVBNQc7GjRsBaGlpobu7u8KlGV/1UEeoj3rWQx2hPuqpOtaOaqhna2vrqLbTqXUREZEqpiAXERGpYgpyERGRKqYgFxERqWIKchERkSqmIBcREaliCnIREZEDlH2+QN+OYkX2rSAXERHZT5kNBR7+t+3c8eebWXVNtiJlqKoJYURERCaDnesG6HIZNv6mFzyYCBQHKjMBqYJcRERklHZ0DdDpenj+7j4ATAzmndpE+9IUTQdXJlIV5CIiInux7cl+Oldk2HR/CPBIHOaf0cTii9IkZ0UrWjYFuYiIyG5seayPzhUZuh/uByCaMCw4q4m2i1I0Tq9sgA9SkIuIiAzhvaf74dAC3/p4CPBY0rDw3CYWXZAiMXVyBPggBbmIiAghwDfdH1rg258aAKAhZVh0QYqF56WIpyfnQC8FuUgdK+VLFFcP4PNgkhBtayCSnJx/rETGiy95nr+nl84VGXauLgAQnxKh7cIUC85poqFpcv+fUJCL1KlSvkThvj68ARMxlPIev6WP2PEJhbnUBV/0bPxdL10uQ8+6EOCJ6RHalqRYcFYTscbq+H+gIBepU8XVAy+GOITfvuQprh4gckSiwqUTGT+lgqfrxm089L3NZJ8Ns7E1tkRYfHGa+ac3EU2YCpdw3yjIReqUz78U4oNMxODzFSqQyDgrDng23JZn1TUZci+EAG86KMriS9LMOyVJpKG6AnyQglykTpkklPL+ZWHuS55Isjr/mInsTrHP88ytOVZdl6G3uwTAlPlxFl2U5JA3JYnEqvvfvIJcpE5F2xrwW/rwJf/iaXXjw+MitaCQL7Hu5hyrr8/Sty0EePP8GO3L0hx1wVy2bttS4RKODQW5SJ2KJCPEjk+82Gs9kjTqtS41YSBXYu0NOVavzDDQE+Y/n9IWo+PSZg4+IYGJGCLR6m6FD6UgF6ljkWREHdukZvT3lFjzsyxrf5plIBsCfNphDXRcmmb2sQmMqZ3wHkpBLiIiVa1vR5HVK7OsuzFHIR8CfMaRcTqWpWl5VbxmA3zQuAe5tbYReD/w18DlzrkrrbXnAn8DNAM3Oec+Nd7lEBGR2tK7tcjq67OsuylHsS8EeMsxcToubWbmEfEKl27iTESL/CAgD/xwyGOPAKcDJWCNtfYK59y2CSiLiIhUufymIl3XZVh/a45SmEmV2ccn6FiWZvph9RPgg4z3E7MQurX2cmCtc+7KIY9NA34PHOGcGxjhNZcBlwE4547t7y9PXh+LUSgUJqDUlVMPdYT6qGc91BHqo56qY2X1PNvPI9/fxKqfb6dUCNm14M1TOPpds5l5WHKf3msy13NQPB4H2Ot1gYpdI7fWpoAfAR8aKcQBnHPLgeXlu767uxuAlpYWBm/XqnqoI9RHPeuhjlAf9VQdKyOzvkDn1Rk23pnHl4AItJ7cSIdN07ygAU+W7u7sPr3nZKzncK2traPariJBbq1tIYT4F5xzd1aiDCIiMrntXDtAl8uw8be94MFEYO4pSdptmvQhk6evdqUXH5qIzm5zgZVAK9BnrT0FmA50AJ+x1gIsd879aLzLIiIik9/2rgE6V/Twwj19AJgYzDs1yeKlaVIHT54Ah8mx+NC4fyLOuQ3AceO9HxERqW7bnuync0WGTfeHAI80wPwzm1h8UZrkrGiFSzeyybD40OT6aiMiInVny6N9dK7I0P2H0KE5mjAsOKuJtiUpGmdMzgAfNBkWH1KQi4jIhPPe0/1waIFvfbw8IilpWHhuE20XpIlPrY6pgifD4kMKchERmTDeezbd10eny7D9qTBgqSFtWHR+ioXnpYinqyPAB02GxYcU5CIiMu58yfP8Pb10rsiwc3UYvx2fEqHtwhQLzmmioam6AnzQZFh8SEEuIiLjxhc9G3/bS5fL0PNMCPDEjAiLl6SY/9YmYo3VGeBDVXrxIQW5iIiMuVLB8+wdebquzpDdWASgsSVC+9I0805vIhqv7YVMJpKCXERExkxxwLPhthDg+U0hwJsOjtK+NM3cU5JEGhTgY01BLiIiB6zY53nm1hyrrsvQ210CIHVIlA6bpvVNSSJRBfh4UZCLiMh+K+RLrLspx+rrs/RtDwHevCBGx7I0c05qxCjAx52CXERE9tlAtsTaG3OsXplhoCesRDZ1cYyOZc0cdEJil0lSZPwoyEVEZNT6e0qs+WmWtT/LMpANAT7tsAY6Lk0z+9gExijAJ5qCXERE9qpve5HVK7Os/XmOYj4E+Iwj4nRcmqblVXEFeAUpyEVEZLd6txZZdV2WdTdlKYWZVGk5Jk7HsjQzj6zc2Gl5iYJcRER2kd9UpOvaDOt/kaMUZlJl9vEJOpalmX5YvLKFk5dRkIuIyIt6nu3nD8u3s+H2PD5MxMbBJzXSYdNMXTxx84fL6CnIRUSEzPoCXddkePbO5/BFIAKtJ4cAb16gAJ/MFOQiInVs59oBulyGjb/tBQ8mCnNPTdJ+SZr0IYqIaqCjJCJSh7Z3DdDlenj+7j4ATAzmndbE8e+fR198Z4VLJ/tCQS4iMo5K+dKLS1yaJBO+xOVw257sp3NFhk33hwCPxGH+GU0svjhNsiVKc0ucvu6KFU/2g4JcRGSclPIlCvf14Q2YiKGU9/gtfcSOT0xomHvv2fJYP10rMnT/IYwhiyYMC85uom1Jisbp0Qkri4w9BbmIyDgprh54McQh/PYlT3H1wISsX+29p/uhfjpX9LD1j2EMWSxpWHhuE20XpIlPrf61wEVBLiIybnyeXeYcNxGDz4/zfr1n0319dK7IsP3pEOANacOi81MsPC9FPK0AryXjHuTW2kbg/cBfA5c756601r4W+BcgCtztnPub8S6HiMhEM0ko5f3LwtyXPJHk+Exn6kue5+7upWtFhp1rwiDw+JQIbRemWHBOEw1NCvBaNBFH9SAgD/xwyGPLgbc7514PHGOtPWkCyiEiMqGibQ0YHwIWwm/jw+NjyRc9z96R586PdPPgl7ezc02BxIwIh7+3mVO+O4v2S9IK8RpmvPcTsiNr7eXAWuAnwK3AmcDXgCOAHzvnvj7Cay4DLgNwzh3b3x86acRiMQqFwoSUu1LqoY5QH/WshzpCfdRzf+pYzBXIP5mjlC0RSUVIvqKJaNPYnAwtFTyrbt7Go9/fzM714e9j6qAGjnrHLNrPnU4sse/hXQ/HEaqjnvF4HGCvp28qdY28FfgB8AngZHZzZsA5t5zQegfw3d1hTERLSwuDt2tVPdQR6qOe9VBHqI967ncd5790szfXD7kDK0dxwLPhlzm6rsmS31QEoOngKO2XpJn7liSRBs/2nq3Qs+/vXQ/HEaqjnq2traPabsKD3Dm3zVqbAd4HPAd8Gfj0RJdDRKTaFPs8z9ySY9V1GXq3lABIHRKlw6ZpfVOSSFRLidajiejsNhdYSWiF91lrTwE+CFwPDAC3OufuG+9yiIhUq0K+xLqbcqy+Pkvf9hDgzQtjdNg0c05qxCjA69q4B7lzbgNw3AhPnTDe+xYRqWYD2RJrb8iy+idZBnpCf6api2N0LGvmoBMSuwxtk/qkceQiIpNMf0+JNT/JsuaGLIVsCPDpr2igfVma2ccmMEYBLi9RkIuITBJ924usXpll7c9zFPMhwGccGefQS9PMPDquAJcRKchFRCqsd0uRVddlWXdzllIYRcasV8dpX9bMzCPilS2cTHoKchGRCsltKrDqmizrf5GjVB7SfNBrE7TbNNMPU4DL6CjIRUQmWPa5Al1XZ9hwex4fhoFz8EmNdCxLM3WMZ32T2qcgFxGZIJn1BTpdhmd/nYcSEIHWNzXScUma5gUKcNk/CnIRkXG2c80AnS7Dc7/rBQ8mCnNPS7L4kjTpVv0ZlgOjf0EiIuNke9cAnVf18MK9fQBEYjDv9CYWX5yi6SD9+ZWxoX9JIiJjbNuT/Tx9VYbND5QDPA7zz2xi8UVpki3RCpdOao2CXERkDHjv2fJoP50rMmx5JIwhizYaFpzdxOILUySmK8BlfCjIRUQOgPeezQ/103lVD9ueGAAg1mRYeG6KtvNTxKdqHXAZXwpyEZH94L3nhd/30bkiw47OEOANacOiC1IsOjdFQ1oBLhNDQS4isg98yfPcXb10rsjQszbM4hKfGqHtwhQLz24i1qQAl4mlIBcRGYVS0bPqlm089N1uMutDgCdmRFi8JMX8tzYRa1SAS2UoyEVE9qBU8Dx7R55OlyH3XJiGLTkryuKlKead1kQ0roVMpLIU5CIiIygOeDb8MkfXNVnym0KANx8SZ9FFSea+JUmkQQEuk4OCXERkiGKv55lbc6y6LkPvlhIA6blR2m2ao5fMY+v2LRUuocjLKchFRIBCvsS6m3Ksuj5L//YQ4M0LY3QsSzPndY2YqCESUytcJh8FuYjUtYFsibU3ZFn9kywDPR6Aqe0NdCxLc9BrE5iIwlsmNwW5iNSl/p4Sa36SZc0NWQrZEODTXxkCfNZrEhijAJfqoCAXkbrSt73I6pVZ1v48RzEfAnzmUXE6lqWZeXRcAS5VR0EuInWhd0uRVddlWXdzllKYCp1Zr0nQYdPMOCJe2cKJHAAFuYjUtNymAquuybL+FzlKYR4XDjohQbtNM/1QBbhUv4oFubW2Efg+0Ao0Al90zl1XqfKISG3JbizQdU2GDbfn8UXAwJzXN9Ju00xta6h08UTGTCVb5CcAfc65N1hrjwX+GVCQi8gB6Vk/QJfL8uyv81ACInDImxtpvyRN83wFuNQe472vyI6ttQngKmALsAD4R+fc74dtcxlwGYBz7tj+/nBhKxaLUSgUJrbAE6we6gj1Uc9armMxVyD/ZI5StkSsuYHEoQmiTZVpH2ztyvPI9zaz9lc7wIOJwuK3TufoP5vFlHmJMdlHLR/LQfVQR6iOesbjcYC99r6sZIt8CpACbiKcXj8eeFmQO+eWA8vLd313dzcALS0tDN6uVfVQR6iPetZqHUv5EoX7+vAGTMSQzCXZvnorseMTRJITt4DI9s5+OldkeOHePgAiMZh7WhPtS1M0HRSjnx66u3vGZF+1eiyHqoc6QnXUs7W1dVTb7TXIrbXTgPcBNznnHh/h+TjwZeBTzrncPpTxfOAR59xya+3PgRuA/9yH14tIBRVXD7wY4hB+exMejxwxNi3gPdn6RD+dV2XY/GA5wOMw/8wmFl+UJtkSHff9i0wWo2mRzwTeAZxprb3fOfeJwSestXOB/wW69jHEAX4O/Km19o5yOT63j68XkQryeXaZ9cxEDD4/jvv0ni2PhQDf8ki41BZtNCw4u4nFF6ZITFeAS/0Z7an1Pzrn3matNeVr28cDbwdeBVzhnLtmX3fsnHsOOGVfXycik4NJQinvXxbmvuSJJMd+QhXvPZsf6qfzqh62PTEAQKzJsOi8FIvOSxGfqrXApX6N+hq5tfYmoAOYD/QDXwQ+7JwrjlPZRGQSi7Y14Lf04UshzH3JY3x4fKx473nh9310rsiwozMEeEOzoe2CFAvPSdGQVoCL7PF/gbX2buBfgdnAUudcOzALOANIAndZay8Y91KKyKQTSUZCx7YZEUgYGmY1jFlHN1/ybPxdnt98rJv7P7+NHZ0DxKdFeMW7mjn1f2bTsaxZIS5Strf/CWcSepLfC9xnrf0H59wOYLtz7p+As4D3W2u/PM7lFJFJKJKM0HBEgvhxCdKvmXLAIV4qejbckefOj3Tz4Je3s3NNgcSMCIe/fwqnfmc27ReniTUpwEWG2uP/COfcTuAC59wngV8Bv7XWHgdcaa1tBj4NfBeYba09aNxLKyIVVcqXGHi8j/77+xh4vI9SvjQ271vwPPOLHHd8cDMPX7GdzPoCyVlRjvzgFE75zmzazk8RbdRiJiIj2Zdx5G8BriVM3mKA1wMPOOeuB64fh7KJyCQyfNx4Ke/xW/oO6HR6sd+z/pc5Vl2bJb8pdLdpmhOl/ZI0c9+cJNKg8BbZmz0GubX2emCxtfY6wnzoQ32svI0FbnDOfXt8iigik8GI48ZLfr/GjRd7PetuzbHq2gx9W0OrPj03SrtN03pykkhUAS4yWntrkf8T8KPy76HzoHvgb4ABYCrwM0BBLlLDxmLceCFXYu1NOVavzNK/PQR488IYHcvSzHldI0YBLrLP9naN/DGgtzyj2/CLYccC/wFkgN+OT/FEZLIwydCbfChf8pjk3l87kCnRuaKH2963iSev7KF/e4mp7Q0c94/TOfnrLbS+IakQF9lPo7lG/ofy778lBHYMeJ1z7vvW2keAjzrnlo5XAUVkctifceP9O0us+WmWNT/LUsiFLwHTX9lAx7I0s16TwBiFt8iBGk2QX2utfQ+wBJjlnPse8BEA59zDwIestTOdc1vGsZwiFVfKl8J14nxonUbbGiZ0cZBKGxw3PvgZRJJmt59B37Yiq1ZmWffzHMXeEOAzj47TsSzNzKPiCnCRMbS3zm4JoAWYAXwS+L61tgfYNGSz3wO3WGsvds6tG7eSilTQePTYrkaRZGSPHdvyW4qsvi7DuptzlMJU6Mx6TYKOZWlmHB6foFKK1Je9tcivBX5H6LH+ZPn3l4HflJ/3QBrIKMSllo1lj+1alNtU4O7/9yydP9tKqbzE80EnJOiwaaYdqgAXGU+7DXJr7Rzg18DtwLuBrwIOWOKce/eQ7f4bLT8qNa4SK31Vg+zGAl3XZNhwex5fBAzMeX0jHcvSTFk0dnOui8ju7facYHl1sn8FDgXeAySAywmTwQAvrkV+KrByXEspUmEH0mO7FvWsH+ChK7bxqw9uZv0v8ngPbWdO403fbOHYf5iuEBeZQHtqkb8auBq4D/i9c+5D5cex1r6F0PktAdzonBuYiMKKVMpErPRVDXauGaBzRYbn7uoFDyYKc09P0r40zYKjD6a7u7vSRRSpO7sNcufcQ0A7gLX2vmFPdwI3EdYlv9RaO6U8L7tITdqXHtu1aPvT/XS6DC/c2wdAJAbzTm9i8cUpmg7al5meRWSs7alFboDDy5PBDP1r5Z1zG4ANwE3W2gHgL4HPjGtJRSpsbz22a9HWP/bTuSLD5gfLAR6HBW9tou2iNMmZ0QqXTkRgz73WZwI/ttZ+Efg7a20r0D/Cdt8m9GxXkIvUAO89Wx4JAb7l0fBfPtpoWHhOE20XpEhMV4CLTCZ7OrXeba09DbgBmAc8CjwCbLTWrilv5oG3ATlr7dxyS11EqpD3ns0P9tG5IsO2J0K3l1jKsOjcFIvOTxGfUh+XEUSqzR4vbjnnNllrzyZMzfrXzrlHR9rOWvsQMJtwul1Eqoj3nhfuDQG+oysEeEOzoe2CFAvPSdGQVoCLTGZ77aVSbpm/AziB0CofyQecc4UxLZmIjCtf8jx3Vy+dKzL0rA3/fePTIixekmLBWU3E6qQjn0i1G1V3U+fcfYRhaC+y1sadc4PXzL9grf2ic27HWBdQRMZWqejZ+Os8XVdnyawPAd44M8Lii9LMP6OJaKPmQRepJnsNcmvtcUDTsIcbgD8D3mmtnQr8BfDZsS+eiIyV0oBnwx15uq7OkHuuCEBydpT2pSnmntZEtEEBLlKNRtMi/z9Ca9wA5xA6v6WBBmvt4cBRwAPOuey+7txa+1pCb/cE8EPn3Hf39T1EZM+K/Z71v8yx6pos+c0hwJvmROmwaQ55c5JITAEuUs1Ge2r9nQDW2kecc++01nYS5l//NNAG/Ne+7tha2wx8HTjfObd5N9tcBlxWLgMtLS2h0LHYi7drVT3UEeqjnpWqY6G3xFMrt/L4D7vJdYdT6FMXJjj6XbNYdOq0MQ9wHcvaUA91hNqq52iC3I/0mHPut9barxNa59/fj32/EYgD/1sO9a855346dAPn3HJg+eA+B6d/bGkF2RUJAAAbnElEQVRpqfmpIOuhjlAf9ZzoOhZyJdbelGP1yiz920sATFkUo31Zmjmva8REimzdvmXM96tjWRvqoY5QHfVsbW0d1XZ7W4/8YsJp75GeOwY4BDjTOVfa1wICU4D7nXN/bq2dRVjX/Kd7eY2I7MZApsSaG7Ks+WmWgZ7w/XtqewMdl6Y56LUJjNEpdJFatKcpWo8idGgbaX2ng4ErgaXAUmvt2v3osd4JvKt8ux/Y52vsIgL9O0us/kmWtTdkKeRCgE9/ZQMdy9LMeo0CXKTW7WkZ00edc+cDIwV0D3CSc+63wHTgrH3dsXPuAeBJa+1dhJb4x/b1PUTqWd+2In/83k5ue+8mulyGQs4z8+g4J35hBid9ZSazj21UiIvUgf29Rp5zzuXKt28mTNN61b7u3Dn38X19jUi9y28psuraDM/ckqNUnslh1msSdCxLM+PweGULJyITbjRBbqy1HyIMP5tZvj21/BvCte63DpsgRkTGWO6FAquuzbL+FzlK5XkUDzohQYdNM+1QBbhIvRpNkH8TmFW+vbx8+z+HPNYArAY6gMfHuoAi9S67sUDX1Rk2/CqPLwIG5ry+kY5laaYsaqh08USkwkYT5Ic653a5fm2t/VvgEefcLcCnxrxkInWu55kBuq7O8Oyve6EEROCQNydptyma5ynARSQYTZC/jWEd0co92j8JnDQehRKpZzvXDNC5IsNzd/WCBxOFuacnaV+aJtU6qjmcRKSOjOoa+dA71toI8B3gk865J8alVCJ1aPvT/XS6DC/c2wdAJAbzTm9i8dIUTbMV4CIysv3ptf4xYIdz7r/HoTwidWfrH/vpXJFh84PlAI/DgremWHxRisaZ0QqXTkQmu9EEedRaO905t618/zDg/eNYJpGa571nyyMhwLc8GgZ7RJOGhWc30XZhisQ0BbiIjM5ogrwBWGOt/TJwK/Ac8C5r7S4bOue0lKnIHnjv2fxAH50uw7YnBgCIpQyLzk2x6PwU8Sm7naNJRGREownyPuDtwLcJq51dC6wdxzKJ1Bxf8rxwbwjwHV0hwBuaDW0XpFh4boqGlAJcRPbPqK6RO+dusdYeSVjl7E3A55xzT49v0USqny961vxyOw9+t5uetWEWl/i0CIuXpFhwVhOx5NgEeClforh6AJ8Hk4RoWwORMXpvEZncRt0V1jmXAS4qn2K/01p7inqti4ysVPRsvDNP19UZMhuKADTOjNB2UZoFZzYRTYzdHOilfInCfX14AyZiKOU9fksfseMTCnOROrDP/8udc/8A/B+w0lqbHvsiiVSv0oDnmVtz3PHBzTz8bzvIbCiSOriBoz40hbd8ZzZt56fGNMSB0BIvhziE396Ex0Wk9u3zOHIA59zfWWvnA/8P2LXXm0idKfZ71v8ix6prs+Q3hxZ405woHTbNq5bOY+v2LeO2b59/KcQHmYjB58dtlyIyiYwmyP9zN4+/D/iWtTZdPu0uUneKvZ51t+RYdV2Gvq0lANLzYnTYNHPe2EgkaojExncpUZOEUt6/LMx9yRNJaglTkXqw1yB3zl2+m8czhN7sInWnkCux9uc5Vq/M0r8jBPiURTHal6WZ87rGXVrI4yna1oDf0ocvhTD3JY/x4XERqX2a91FkHwxkSqy5Icuan2YZ6AmTHk7taODQS9PMPj6BMRPfCo4kI8SOT7zYaz2SNOq1LlJHFOQio9C/o8Tqn2ZZe0OWQi4E+IzDG+hY1kzLq+MVCfChIskIkSMSFS2DiFSGglxkD3q3FVl9fZZ1N+Uo9oYAn3l0nI5L08w8svIBLiKiIBcZQb67yKprMzxza45SmAqdWccm6FiWZsYr45UtnIjIEApykSFyzxfouibL+tty+DARGwedmKDDppnWoQAXkclHQS4CZDYW6HIZnv1VHl8CDMx5QyMdNs2URer9LSKTl4Jc6lrPugE6r86w8Te9UAITgTlHx2k7LcnUExvV81tEJr2KBrm19iPAvwPtzrm1lSyL1JcdqwfoXJHh+bt7wYOJQuvRcRaekCA1M4YveQr3ab5yEZn8Khbk1tpFwHnAXZUqg9SfbU/307Uiwwu/7wMg0gDzTm9iwdENJErm5fOVlzzF1QMa1iUik1pFgtxaa4D/AD4KLN/DdpcBlwE452hpaQEgFou9eLtW1UMdYeLq+cIfsvzhe5vYeG+YTTiaMBy2ZAZH/sksmmY1sPM32yn1lnZ5XSQSYUrLtAPat45l7VAda0ct1bNSLfLLgNudc09bu/s1V5xzy3kp6H13dzcALS0tDN6uVfVQRxjfenrv2fJIP09flWHrY2EMWTRpWHhOE20XpEhMi5JjB7luGCj1UcqUdp2vPB6hv7twQOXQsawdqmPtqIZ6tra2jmq7SgX5ecAUa+35wDHAVdbaP3POPVWh8kgN8d6z+YE+Oldk2PZkWMozljIsOi/FovNTxJt3veat+cpFpFpVJMidc+cO3rbW3gG8S53d5ED5kueFe/vodBl2dIUAb2g2tF2QYuG5KRpSu++0pvnKRaRaVXz4mXPuzZUug1Q3X/Q8d1cvnS5Dz9pwGjwxLULbkhQLzmoiNsowrvb5ykv50otfREwSfRERqRMVD3KR/VUqejbemafr6gyZDUUAGmdGWHxRmvlnNhFN1M886KV8icJ9fXgTetyX8h6/RcPnROqBglyqTmnAs+H2PF3XZMg9HwI8OTtK+9IUc09rItpQPwE+qLh64MUQBw2fE6knCnKpGsV+z/pf5Fh1bZb85hDgqdYo7ZekOeTNSSKx+gvwQT7Py3rcQznM8xUqkIhMGAW57NF4X3ct5goMPN63x/cv9JZ45uYcq67P0rc1jPVOz4vRYdPMeWMjkWj9Bvggk4RS3u86fC6pz0ak1inIZbfG+7prKV9i531hIpaR3r+QK7H2xhyrf5Klf0cI8CmLYrQvSzPndY27tEDrmYbPidQvBbns1nhfdy2uHqAh0rjL+/c+3seGzgJrfpplIOMBmNrRwKGXppl9fAJjFODDaficSP1SkMuISvkShTUFyHuIgZkeIdIQGdPrrj4PZshp8f5ciXX39bH+wT4KfSHAZxzeQMeyZlpeHVeA70W1D58Tkf2jIJddvHhKvb8EA4SfXJHSISF4x+q6q0mC7/X0ZUqs/X0vGx7qoxjmcaHlVXE6lqWZeZSCSURkTxTksovBU+rRGTFKGwt4wBtga4lIS3TMrrv2T43w9HXbWfdgjlJ5OvOWthgd72lm5qsax2QfIiK1TkEuu3hxKFMEIq0xStuKUADiZkw6uuVeKNB1dZb1t+Xw5QCf/YoG2k5LMuPkpK7riojsAwW57GLoUCbTYIjOjoWhTDMiBxSymY0FulyGZ3+Vx5cAAwtPm8r8CxqYslC9q0VE9kddBrnmpN6zsR7K1LNugE6XYeNve6EEJgKHvCVJxyVpFr764Em/lKCIyGRWd0GuOamDPX2ZGauhTDtWDdC5IsPzd/cCYKIw94wk7UvTpObU3T89EZFxUXd/TTUn9ei+zBzIUKZtT/fTeVWGTff1hfdqgHmnN9F+cZrk7OiY1UNEROowyDUn9fh9mdnyeD+dV/XQ/XA/AJE4LDgrxeIlKRpnKsBFRMZD3QW55qQe2y8z3nu2PNLP01dl2PpYCPBo0rDw7CbaLkyRmKYAFxEZT3UX5JqTemy+zHjv2fRAH51XZdj+VJjFJZYyLDovxaLzU8Sb66e/gYhIJdVdkGtO6gP7MuNLnufv7aNrRQ87VoVB4A3NhrYL0yw8p4mGVP18jiIik0HdBTloTur9+TLji57n7uqlc0WGnnUhwBPTIrQtSbHgrCZidfRFSERkMqnLIJfRf5kpFT0b78zT6TJkny0C0DgzwuKL08w/o4loon76FoiITEYKchlRacCz/vY8q67JkHs+BHhydpT2S9LMPTVJtEEBLiIyGSjI5WWK/Z71v8jRdU2G3u4SAKnWEOCHvDlJJKYAFxGZTBTkAkCht8QzN+dYdV2Wvm0hwJvnx2i3aVrf0PiydcNFRGTyqFiQW2ubgOXAIUAT8E7n3FOVKk+9GsiVWHdjjtUrs/TvDAE+pS1Gx7I0B5/YuMt4cxERmVyM975iO7fWLnbOrbLWXgYc4Zz72LDnLwMuA3DOHdvfHyYcicViFAqFCS/vRBrvOvbtLPLE1d38ccUW+nvCNfCWw5O86t2zmfv6ZoyZmADXsawd9VBP1bF2VEM94/E4wF7/GFc0yAdZay8Htjnnvr6HzfzGjRsBaGlpqfkVs8arjn07iqz5SZa1N+Qo5MOxn3FEnI5laVqOiU9YgA/Ssawd9VBP1bF2VEM9W1tbYRRBXvFr5NbaC4FXAxdXuiy1rHdrkdXXZ1l3U45iXwjwlmNCgM88sn7H1IuIVLuKBrm19r3AKcAlzrnJfY6jSuU3F1l1XYZnbs1RClcmmH1cgo5laaa/Ij4hZdif9d+1ZryIyOhUsrPbcYTObr8DbrXW9jvnzqhUeWpN7vkCXddkWX9bDl/+inTw6xJ02Gamtk/cvPL7s/671owXERm9igW5c+5+QEtjjbHMswW6XIZn78jjS4CB1pMbabdppiyY+IVh9mfJVK0ZLyIyehW/Ri5jo2fdAJ0uw8bf9kIJTATmnpKk/ZI06bmVO8z7s2Sq1owXERk9BXmV27FqgM4VGZ6/uxcAE4N5pydZvDRN6uDKH979WTJVa8aLiIxe5f/Sy37Z9lQ/nVdl2HR/HwCRBph3RhPtF6VJzp48Vyz2Z8lUrRkvIjJ6CvIqs+WxPjpXZOh+OHRBjyYMC85qom1JisYZkyfAB+3PkqnDX0PUAxEKjw+oB7uIyDAK8irgvaf7D6EFvvXxcoAnDQvPaaLtghSJaZMvwIfan/XfB1/zUg/2EhSh8EyBwmMDRA6N0vAK9WIXEVGQT2LeezbdH1rg258aAKAhZVh4fopF56WIN9d+iA32YKcIpY1hHJ3HU1pTpNCjIWkiIgryPajUpCS+5Fl3xw4e/J9udqwK4dXQbGi7MM3Cc5poSFVXcB3I5zjYg73YXcDDS1PIFsEbNCRNROqegnw3KjEpiS96nrurl84VGXrWhQBPTIvQtiTFgrOaiO3Hfis9Q9qBfo6DPdgpvBTi3ntMzGhImogICvLdmshJSUpFz8Y783S6DNlnw0pkTbMbWLQkyfzTm4gm9m/Y1WSYIe1AP8cXe7BHPYSrCxgPZnpEQ9JERFCQ79ZETEpSGvCsvz3Pqqsz5F4IAZ6cHaX9kjTH2Lls27l1968dRUt7MsyQdqCf42APdp6MUHy6HxoMZkYEEzUakiYigoJ8t8ZzUpJiv2f9L3J0XZOht7sEQKo1SrtNc8ibkkRihmh89y3m0ba0J8MMaWPxOUaSEeKvTlB6RYMWUhERGUZBvhvjMSlJobfEMzfnWHVdlr5tIcCb58doX5am9fWNmOjowm20Le3JMEPaWH6O+zOMTUSk1inId2N/JjLZnYFciXU35li9Mkv/zhDgU9pidCxLc/CJjbu0mvdmtC3tyTBD2lh+jiIisisF+R4caAuwP1Ni7c+yrPlploGMB2DaYQ10LEsz+7jES0Op9tFoW9qTJUTVkhYRGT8K8nHQt6PImp9kWXtDjkI+BPiMI+J0LEvTckx8vwN80L60tBWiIiK1TUE+hnq3Fll9fZZ1N+Uo9oUAbzkmBPjMI8cuTMeqpV3pMeYiInLgFORjIL+5yKrrMjxzS45Seazz7OMSdCxLM/0V8XHZ555a2qMJ6MkwxlxERA6cgvwAZJ8vsOqaDOtvy+PDRGwc/LoEHbaZqe2VGd882oCeDGPMRUTkwNVdkI/F6eTMhgJdV2d49o48vgQYaD25kfZL0kxZWNkJSkYb0JNhjLmIiBy4ugryAz2dvHPdAF0uw8bf9IIHE4G5pyRpvyRNeu7k+ChHG9CTYYy5iIgcuMmRPhNkf08n7+gaoNP18PzdfeF1MZh3apLFS9OkDp5cH+FoA3oyjDEXEZEDN7lSaJzt6+nkbU/207kiw6b7Q4BHGmD+mU0sXpImOTs63sXdL6MN6MkyxlxERA5MRYLcWhsFvgEcCUSBDzrnHh3v/Y62tbrlsT46V2TofrgfgGjCsOCsJtqWpGicMTkDfNC+BLTGmIuIVL9KtciXAhHn3JustW8ArgDOGO+d7qm16r2n++HQAt/6eAjwWNKw8NwmFl2QIjF1cgf4UApoEZH6Ybz3E75Ta+3XgVspt8aBw51zC0bY7jLgMgDn3LH9/eWAjcUoFAr7te9irkD+yRylbIlIKkLjYUmeezjPI1duYvNj4Rx7vDnCK20Lh18yk8TUynzXOZA6VpN6qGc91BHqo56qY+2ohnrG43GAvfZAruQ18k8BNwHnA0+OtIFzbjmwvHzXd3d3A9DS0sLg7f0yP5xSf/6eXjqveJadq8PBjE+J0HZhigXnNNHQFKFnYDs9B7CbA3HAdawS9VDPeqgj1Ec9VcfaUQ31bG1tHdV2lQryB4AG59xnrbWvA56YqB37omfj73rpchl61oUAT0yL0HZRigVnNRFrVGcvERGpHpUK8h8DZ1hr7wQKwAcmasebHujjoa9tB6CxJcLii9PMP72JaELjp0VEpPpUJMidcwPA2yux79nHJZh9fIKDTmhk7ilJog0KcBERqV51NY4cwrjx1356RqWLISIiMiZ0QVhERKSKKchFRESqmIJcRESkiinIRUREqpiCXEREpIopyEVERKqYglxERKSKVWTRlP1UNQUVEREZI3udtayaWuRm8Mda+8DQ+7X4Uw91rJd61kMd66WeqmPt/FRRPfeqmoJcREREhlGQi4iIVLFqDfLle9+k6tVDHaE+6lkPdYT6qKfqWDtqpp7V1NlNREREhqnWFrmIiIigIBcREalqCnIREZEqFqt0AUbLWhsFvgEcCUSBDzrnHq1sqcaWtbYXuKd89+vAs8C/EOp7t3PubypVtv1lrW0E3g/8NXC5c+5Ka+1rGVavaj6+u6nj5cClwPPADufcBdVcRwBrbROhg9AhQBPwTmAqtXUsR6rj26ixYwlgrb0FiAONwGVAkho6ljBiHS+mBo9lNbXIlwIR59ybgH8ArqhwecbD8865N5d/rif8QXm7c+71wDHW2pMqXL79cRCQB3445LGR6lXNx3ekOgJ8uXwsLyjfr+Y64pzLAf/snHsL8F3gQ9TYsdxNHaHGjiWAc+7Mcj1vAk6ixo4ljFhHqMFjWTW91q21XwdupfyNCTjcObegsqUaW9baLmAj0A18AfgWcCbwNeAI4MfOua9XroT7r9xCXQv8hHAcX1YvoI0qP76DdSy3yD8OvJ0Q8Fc7575RS/+Gy3X1wLnU4LGEF+u4jTC7Vs0dS2vtOcBngGnAm4CV1NixHFbHkwFLDR7LamqRA3wKOAY4HyhUuCxjzjnX7pw7GfgB8HmgtXz7G+za2qtmu6tXzRxf59y/O+eOA04H3m6tPaz8VNXX0Vp7IfBq4L+o0WM5pI7/WavH0jl3Y7lef0FoNNTcsRxex1o9ltUU5A8ADznnPgscBzxR4fKMJwOsATLA+4BHgLOBuypZqLHgnNvGyPWq1eNbLP/uoQbqaK19L3AJcIlzbjM1eCyH1XHoH/WaOpZDbAdy1OCxHGI7oRU+qKaOZTWdWm8AvgfMI3xj+oBzrrOypRo71trZwLVAH7AJ+DChRfAlYAC4tfwPrapYa+cSTtm1Eur2G+BKhtWrmo/vbur4BHAG4fTzt5xzrprrCGCtPQ64F/gdUAL6gS9TW8dypDr+ito7losI/w97gZ3Ax4HDqK1jOVId30mNHUuooiAXERGRXVXTqXUREREZRkEuIiJSxRTkIiIiVUxBLiIiUsUU5CIiIlVMQS5SQ6y1s6y1kfLt15R/t1prTWVLJiLjRUEuUiPKcxHcTpheE+Bqa20ceAdwo7U2OWTbI621K4fcn2qtvYe9sNY2WGt/Y619vDz+dnfbXWitPW2/KyMio6Zx5CI1wFrbQlgY4t+ccz8qP3YV8DPgacKUlOuBywnzhzcSVvhaVX6LlcAJwF8Oedst5Zn4BvcRAb5DmLd6O6Eh8F7nXGnINh8gzGXdTJhg44/AFqCBMGHOQsKERzngHufcB8boIxCpW2qRi1S58rKwvwPmA5uttQ9aax8gBPV3gfcQlqicA3ypPNf0pcCdzrnjyvdTQJowx//ngesIS7MO7iMJ/B9h9q+3A38OzAZ+VF76c9Acwhzs/+ucm+GcewMwl7BYxWnAOsL0n6cBHxuHj0Ok7qhFLlLlrLWXEaYVvQJ4l3Nuw5DnvkWYNvbvCCH8ENBJ+BIfJ0xfeTtwCnA8cI1z7kxr7S+Bv3TOPWqtPZKwiMZq4E+cc/nyezcQWuivAy5zzt05ZJ/nEla7g7AgxfsIUw4Pznd9EPAJ59wPxvwDEakzsUoXQEQOjHNuOYC19g+E0B7qQ4TT6YPPfQx4ctg2/cAGwqIZ7dbaKPBAOcT/Dvgo4RT92cAT1to04UvA1vLrfwwst9ZeRzgzsBA4q/x+PyEswgHwp8DBhC8E5x9ovUUkUJCL1I4jgV9Za0d6Lg58GrgRWDLsuaecc1+01k4jnJ6PO+f+vvzcCuCbzrnc4MbW2vcBRzrnPj7ksU8RTs2fTFicYnm5PJ93zm0dUqbTgNsOqJYi8jIKcpHacZhzbiGAtXY+0OGcu618/0pCJ7OvEjq1vYy19l3AVMLKV8cDvy4/dRLwtWFfDpqABmvt0mFv87fAM8A/AJsJre73W2uvB7LlbRKE0/kiMkbU2U2kNh0JXDj8Qefch4H/AN4NdBNa6ac65x4jnA7/IuEU+OD2P3bOzQUWA28s3/474Lvl229xzs0t//yY0KJ/JaGj3b8Qrsv/PS+ZTejFLiJjREEuUpuOBxqttYmhD1prjyH0Oh8oPzQL+Jy1diZwJvAV4Jhyi36oLwJfG/ZeU4FfW2tfFvzAgHPuROfciYTr8RHC9XRPaOEXgB5eaqWLyAFQkIvUjkuttTOttZ8E3kYYq/2YtfbbwDWEnu0fIFzDbgaKhCFldwDfJwxNGyC0pK+01jYCWGvPBZYROs69yDm3g9AB7hvW2vfspWwfA94FPAVcDbzgnLv+QCssIgpykapXnm3tu4Tx23cRrkOf6Jz7GOEU++2E69e/A/4J+Cyh09uN5clczgMed85dB+Ccu5owUcyXrbUnAFcCS51zm8q7nEXokY5z7iFCmF9hrT2j/HzUWntPeaa4cwgzzd1C6Ax3DvAGwpeOz43TRyJSV9TZTaTKOecGrLX/DWxyzj0z7Lk+Qs/zFdbaRc65zcDFw97ik86554Y99ueETm2NgHXO3WOt/RfC9fNeQvgP7uNea+2pwBPlh85xzg0OTcNaOw9Y5Zx7uPxQHyHIF+5/rUVkkCaEERERqWI6tS4iIlLFFOQiIiJVTEEuIiJSxRTkIiIiVUxBLiIiUsX+P2NCyHmHNZ8eAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 画出预测散点图\n",
    "ax.plot(xs_pred, ys_pred, '-', color='darkorchid', linewidth=2)\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>        <td>不良贷款\n",
       "(亿元）</td>    <th>  R-squared:         </th> <td>   0.712</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.699</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   56.75</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Thu, 20 Dec 2018</td> <th>  Prob (F-statistic):</th> <td>1.18e-07</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>22:13:25</td>     <th>  Log-Likelihood:    </th> <td> -51.508</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    25</td>      <th>  AIC:               </th> <td>   107.0</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    23</td>      <th>  BIC:               </th> <td>   109.5</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "       <td></td>          <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>       <td>   -0.8295</td> <td>    0.723</td> <td>   -1.147</td> <td> 0.263</td> <td>   -2.325</td> <td>    0.666</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>各项贷款余额\n",
       "(亿元)</th> <td>    0.0379</td> <td>    0.005</td> <td>    7.534</td> <td> 0.000</td> <td>    0.027</td> <td>    0.048</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>14.277</td> <th>  Durbin-Watson:     </th> <td>   2.464</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.001</td> <th>  Jarque-Bera (JB):  </th> <td>  14.611</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 1.382</td> <th>  Prob(JB):          </th> <td>0.000672</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 5.527</td> <th>  Cond. No.          </th> <td>    262.</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:              不良贷款\n",
       "(亿元）   R-squared:                       0.712\n",
       "Model:                            OLS   Adj. R-squared:                  0.699\n",
       "Method:                 Least Squares   F-statistic:                     56.75\n",
       "Date:                Thu, 20 Dec 2018   Prob (F-statistic):           1.18e-07\n",
       "Time:                        22:13:25   Log-Likelihood:                -51.508\n",
       "No. Observations:                  25   AIC:                             107.0\n",
       "Df Residuals:                      23   BIC:                             109.5\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "===============================================================================\n",
       "                  coef    std err          t      P>|t|      [0.025      0.975]\n",
       "-------------------------------------------------------------------------------\n",
       "const          -0.8295      0.723     -1.147      0.263      -2.325       0.666\n",
       "各项贷款余额\n",
       "(亿元)     0.0379      0.005      7.534      0.000       0.027       0.048\n",
       "==============================================================================\n",
       "Omnibus:                       14.277   Durbin-Watson:                   2.464\n",
       "Prob(Omnibus):                  0.001   Jarque-Bera (JB):               14.611\n",
       "Skew:                           1.382   Prob(JB):                     0.000672\n",
       "Kurtosis:                       5.527   Cond. No.                         262.\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 显示统计检验结果\n",
    "# fitted_ols.summary2(float_format=\"%.4f\")\n",
    "fitted_ols.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "--------------------\n",
    "\n",
    "判定系数: $R^2 = \\dfrac{SSR}{SST} = 0.712$, 拟合优度的度量, 在不良贷款的取值变动中有71.2%是由贷款余额所决定的.\n",
    "\n",
    "多重判定系数(调整后): 解释变量变多, 势必导致误差SSE变小, 进而$R^2$变大(高估了), 修正$R_\\alpha^2 = 1 - (1 - R^2)(\\dfrac{n-1}{n - k - 1}) = 0.699$\n",
    "\n",
    "$\\beta_0 = -0.8295$, P-Value($\\beta_0$) = 0.2631, 置信区间[-2.3252, 0.6662]\n",
    "\n",
    "$\\beta_1 = 0.0379$, P-Value($\\beta_1$) = 0.0000, 置信区间[0.0275, 0.0483]\n",
    "\n",
    "### 1. 线性关系的F检验\n",
    "\n",
    "$H_0: \\beta_1 = 0$\n",
    "\n",
    "$F = \\dfrac{SSR/1}{SSE/(n-2)} = \\dfrac{MSR}{MSE} ~ F(1, n-2)$\n",
    "\n",
    "查表$F_\\alpha(1, 23) = 4.28 < F-statistic = 56.75$, Sig.F(P-Value) = 1.18e-07 < $\\alpha = 0.05$, 决绝$H_0$\n",
    "\n",
    "\n",
    "### 2. 回归系数的t检验\n",
    "\n",
    "$H_0: \\beta_1 = 0$\n",
    "\n",
    "$t = 7.534 > t_0.025 = 2.0687$， P-Value = 0.000 < $\\alpha$ = 0.05, 决绝$H_0$\n",
    "\n",
    "**所有统计量检验都是在$H_0$成立的假设条件下, 计算概率的**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. 区间估计"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 置信区间估计 (y平均值的置信区间估计)\n",
    "\n",
    "$$\n",
    "\\hat{y_0} \\pm t_{\\alpha/2}s_e \\sqrt{\\dfrac{1}{n} + \\dfrac{(x_0 - \\bar{x})^2}{\\sum_{i = 1}^{n}(x_i - \\bar{x})^2}}\n",
    "$$\n",
    "\n",
    "There is a 95 percent probability that the true regression line for the population lies within the confidence interval for our estimate of the regression line calculated from the sample data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.0686576104190406 1.979947532717681 154933.5744\n"
     ]
    }
   ],
   "source": [
    "n = len(xs_loans_balance) # sample count\n",
    "dof = n - fitted_ols.df_model - 1\n",
    "\n",
    "y_hat = fitted_ols.predict(xs_loans_balance_with_intercept)\n",
    "y_err = ys_loans_bad - y_hat\n",
    "x_bar = xs_loans_balance.mean() # xs_loans_balance_with_intercept.T.iloc[1].mean() \n",
    "s_var = np.sum(np.power(y_err, 2))\n",
    "s_err = np.sqrt(s_var/dof)\n",
    "\n",
    "t_half_alpha = stats.t.ppf(1-0.025, df = dof)\n",
    "\n",
    "# TODO: x_var = np.sum(np.power(xs_pred - x_bar, 2))\n",
    "x_var = np.sum(np.power(xs_loans_balance - x_bar, 2))\n",
    "\n",
    "# 计算每个x对应的y平均值的置信区间\n",
    "\n",
    "dalta_conf = t_half_alpha * s_err * np.sqrt(1.0/n + np.power(xs_pred - x_bar, 2)/x_var)\n",
    "\n",
    "lower_95_percent = ys_pred - dalta_conf \n",
    "upper_95_percent = ys_pred + dalta_conf\n",
    "\n",
    "print(t_half_alpha, s_err, x_var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAEkCAYAAAAhCDhSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmcW1l54P3fke6VSiWVanW57V4wTe/tpjcXW8KQkAlkG5IQckPekHWghy0DCWRmYEjSJGEJSybwJplA3gy8CQn5HEJgZphheRPesEPK7r3pbrrbmF7crnKt2qW7nPnjXqnkcrlWqbQ938/HH5ekK+kcWeXn3uc85xxljEEIIYQQgyPW6QYIIYQQYn9J8BdCCCEGjAR/IYQQYsBI8BdCCCEGjAR/IYQQYsBI8BdCCCEGjAR/IbqMUuqZSim1jeMOKKVi0c+3RH8f3s5zO0EpdTD621JK3RD9fHFnWyXEYJLgL0T3eQfwM5sdoJSaBr4IXB/d9QmlVAL4ReB/KaVS644/qpT6dNPtUaXUN1vb7E3bexXwJaXUIcL/d+pt+c9KqY+sP2FRSv2EUuqPm27foJTS23ifrFLqIaXUF7Y47tVKqaM774kQ/UHJIj9CtI5S6u3A6zZ4KAm4QNB038eNMb+ulHo98Pqm+8eAMlBtuu+XjDH/Er3HFPBZ4L8YY/42uu/vgP8JfAf4YeBxYIHwJOImYAi4GHg0er1PA88GfqPpPRaNMcs77fNWlFKXA/9f1IevRYH+68Bbo37+MPAY8Eh030FgFMhG/QD438A08L6mlz5jjCk0vc8Q8CnCz+Aa4JvGmN9d15bbgX8NTEXv8ShwP+FnAXAVcBLwgH8wxvze3j8BIbqP1ekGCNFPomCzPuD8CKCBq4wxCxs8bQr4E2PMn2z0mkqpjwHD0c/PAv6a8AThrFLqDsAAJeAvgY8QBv1DQNIY88roeUeBPzDG/FR0+11ABviD6G2uAz4GvGd3Pd+YUurHgA8C48BQU3t94O+BjwOLwKXAsjHmx6Pn/QTwr40xb4xufxwYaWrvDPAWws8VpdQ48AlgmfCEJg38/0qpEeC3jDF+9LwjwB8Bw8aYjymlRglPGF4AJIBvAT8YtXGllZ+FEN1Egr8QbaSUygAfAv7wAoF/p24CHOD9wAPGmFua3uvPgRphgLwaeJdSag54L2GqPaGUepBwuOCFhAH0740xL1ZK/SNhNqHVrgd+CPi8MeafgHptggL+F+EV//uBHwderJQ6DLwRiAPx6MTp01G/fxB4rzHmF5VS3wb+KXqt5xOeEP1v4NeNMQGQV0r9IOEJwdeVUrcZY+42xvyKUupzwBVKqVdHbfy3UVvSwFL0OVwGvBT4ahs+EyE6ToK/EG0SBbi/Ikxjf36Lw9/aFIzWuwT4fwCMMR+OXvtuzh0WAHgtcDtQf+z/Jky3v37dcTXgCaBAGATjwAljzL1b92pnjDHvjdr7L+vuN0qplxJehd9FOCzyWuAEYfq/mQV8l3Do4pLoKv9vjDGLSqkPAj9KeCLwY8CPKaXGCIdXctHzPwJ8Rin1DsKTimXCE5Jp4E+i14bw5OJHCE+OfqclH4AQXUqCvxDt8y7gWcAngZ9QSr0KeAj4QFMauu6dW6T91ztKmNbe6CkJ4HeMMVoplQYuWvf4Q8aYd0ZB8jIgYYz5j1t1Rin1k8A7gc8ZY94UBeHfMMZsJ1DeopS67wKPjQI/Y4z5l2hGwPr2njDGfE4pdRPh0MkyYVEkhEMKv2mM8Zra+QfAgjGmuWDwDwmzHz8RPfd/EA4B/Joxptb0Ob6M8ORJiL4mwV+INlBKvQf4ZcKx5FdEd78L+DDwcqXUK40x9zQ95XeVUm8kHNsfAeaaHrsPsNe9xdXGmCPRe10GXBml1VFKfRSYj457D2vFbM3t+xXCoOsSpv+/vI1uvZuwWO73lFJvIUyJb7diPtPU3qNAyhgzG93+56b2/jnh+P/69r6YcAhhXCl1mTHmseihnwbesO4kKAsESqk3r3uZlwN54O3AceDXgd9USv1A0zFJoLLNPgnRs2SqnxAtpJQaiSrvfxr4IWPMg/XHjDGnjDEvAj4KfEUp1Xy1/XZjzBWE1ebfBZzo9q8SBv5fjqrZN3IU+KmNHjDGvI4w/f+rhIWAvxO16z7CdPk7gV/YZve+YIx5kjA9fyNhynw31fD/KvqzUXt/mvAE4G2EAfoPgecYYxY3aq8x5r3GmEui9twc/fxB4Pbo5x8wxlwS/fkq4RDKtYQnEu8hzDK8takJ04QFiEL0NQn+QrTWnxNm1G6NAux5jDF/BjwPuGeDxwqEUwU/G1XGvx74j8aYVxhjLnRFOkNYSZ9c/0CUKn8F4RU+wAHg95VSk8CLCYPrTVH2YFPGmDdEf1eNMS83xtxojLlrq+ddoL2ZqNZgfXtfAjyn6a7rgH+vlLqScAbDHxFmTrLrnvph4JyhC6XUFcA3lFIvbOrDXwBLxpjnGGOeQ3iV7xJmHkaBKwhPtpY5v6ZCiL4hwV+I1nqNMeZlxpjcZgcZY+43xtSr621gRCn1b5RSf0pYdPYjhOn/ZwBHo/H5Zi9XSk0qpd4K/DzhVL/7lFIfIpxC963ouFcTVs+PEE6v+2vgnwkLEd9ljHEJ585/dJPMQiu8XCl1sVLq/YRX3ZcBD0S3/4CwABHCoZL/3NTe9xPOu/8Y8BZjTCn6+c/U2uqGrwWuBH67+Q2NMY8QLnr0qWjYYDM/R/g5fA34EmFdw+zeuixEFzPGyB/5I3/a+IcwuN2+yePHCa80Pw78X8B4dL8F/BphQKoSTtmzCefz30FYPPh2YCw6PkkYxL4UPX4guv+ThFe2r41ufwR4z7o2/AXwx23o+zjwd4R1C/cAv0k43g/hOgOvJKzw/yJh4SHAN6P2viS6/UXgdU2vGQc+R3hS81PAU8Azmh7/EPDKptsvIaz8vyG6fSp6j29G978w+vn3AUU4NPB54N91+rsjf+RPu/7ICn9CtFlUfe4ZY26/wOM3AA+bC6f1iRarsYwxy0qpY8C8WSt62+j4pxtjvnuBxw4ZY55ad1+McOGbwkbP2Qul1PcBjxpjzmxyzOXGmJMXeGyj9g4RTue7BJgwxhyPlv99PuFc/R9qfj+l1HOBbxljAqXUhDFmqemx64CaCTMF9ftiwMXGmPoKg0L0FQn+QgghxICRMX8hhBBiwEjwF0IIIQaMBH8hhBBiwEjwF0IIIQaMBH8hhBBiwEjwF0IIIQaMBH8hhBBiwEjwF0IIIQaMBH8hhBBiwEjwF0IIIQaMBH8hhBBiwEjwF0IIIQaMBH8hhBBiwEjwF0IIIQaM1ekGtJHsVSyEEGIQqa0O6Ofgz+nTpxs/T01NsbCw0MHWtN8g9BEGo5/Sx/4xCP0chD5Cb/Tz8OHD2zpO0v5CCCHEgJHgL4QQQgwYCf5CCCHEgJHgL4QQQgwYCf5CCCHEgJHgL4QQQgwYCf5CCCHEgJHgL4QQQgwYCf5CCCHEPqtWq5w6dYpyudyR9+/rFf6EEEKIblIul5mbm6NYLBIEAVNTUx1phwR/IYQQos0KhQLz8/OUy2Usy8K2bWq1WsfaI8FfCCGEaANjDPl8nvn5earVaiPodwMJ/kIIIUQLGWNYXl5mcXER13W7KujXSfAXQgghWiAIAhYXF1laWsL3fSzLwrK6M8x2Z6uEEEKIHuH7PmfPnmVlZQVjDPF4vGuDfl13t04IIYToUp7nMTc3Ry6XAyAej3e4Rdu3b8HfcZwh4FXAm4DbtdYfdRznduDlwBlgVWv9k5s8/1nA+4A48A2t9Zvb32ohhBDiXLVajTNnzlAoFIjFYj0V9Ov288r/IFAG/mbd/e/WWn90G8//MPASrfVjjuP8o+M4z9Naf735AMdxbgNuA9BanzN/0rKsjs2n3C+D0EcYjH5KH/vHIPRzEPoIYdBfXFwkn89jWRZjY2N7er1qtcrExATZbLZFLdw+ZYzZ1zeMrvZPRVf+bwReQXhS8Amt9Qcv8Jxx4AvAi4H3AtcDH9daf2CTtzKnT59u3JiammJhYaE1nehSg9BHGIx+Sh/7xyD0s9/7WCwWmZubIxaLUa1WUUq15HVrtRpHjhwhk8m05PUADh8+DLBlAzu6vK/W+o+11seAHwZe4TjO1Zscfhj4GPBBzs8eCCGEEC1jjCGXy/HII49w6tQpXNclkUi0LPB3Wres7e9Hf+cBHMeJO47zWcdxfghAa70MFIBXAvcAPwZ8faMXEkIIIXbLGMPS0hKPPPIITzzxBEEQYNt23wT9uv0s+LsE+DThFXzVcZwXAg8ALwIM8D6tdT1PbwHXAgeaXuI1wKcAF/iC1np2v9ouhBCiv/XSHP1W2Leeaa2fAI5t8NC7Nji2ChxZd98XgWe3pXFCCCEGUi/O0W+F/u+hEEIIsU4vz9FvBQn+QgghBka1WmVubq6n5+i3ggR/IYQQfa9cLjM3N0exWOz78fztGOzeCyGE6GuFQoH5+XnK5XJX7q7XKRL8hRBC9JX6HP2zZ89Sq9WIx+MS9NeR4C+EEKIv1OfoLy4u4nmepPc3IZ+KEEKInhYEAQsLCywvLw/EHP1WkE9HCCFET/J9n/n5eVZXVwmCQIL+DsinJIQQoqe4rtuYox+LxRp/xPZJ8BdCCNETqtUqZ86coVgsEovF5Cp/D+STE0II0dVKpRJzc3OUSiVJ7beIfIJCCCG6Uj6fZ35+nkqlInP0W0yCvxBCiK5hjGFlZYWFhQVc15Wg3yYS/IUQQnRcEASNOfoyXa/95JMVQgjRMb7vN+boD9KWurUFxfK3Uky92Cdz3f6/f/9/wkIIIbrOIG6pG1Qgf7dF7rhN+WQYfp9KuFwkwV8IIUQ/q9VqzM3Nkc/nB2JLXRNA+dE4q7M2hXstTE0BoGzD8PU1Djw71ZF2SfAXQgjRdoO2pW5tQZGbtcmdsPGW1xYgSl3ukT3mkrnRw4/VGDsy1ZH29fenL4QQoqMGaUtdvwz5u21ysxaVU2vh1ZoIyB5zyR5zSUyateNrnWhl1KbOvbUQQoh+NEhb6poASg/HydXT+l6U1k8YRm4Mr/JTl/uoLlt9WIK/EGJHgnKAf9LFlEGlIH65TSzVZf+ziY4YpC11a3MxVo9b5E/YeKtNaf0rwoA/8kyPWLKDDdzCvv2rOI4zBLwKeBNwu9b6o47j/ATwZmAE+KzW+m2bPP924OXAGWBVa/2T7W+1EKJZUA7wZqsYBSqmCMoGs1jFmknKCcAAq++ut7y8TBAEfTtdzyvCytdtcsdtKt9bK1S0J9fS+vaE2eQVusd+/uscBMrA3zTddw/ww0AAfNdxnPdrrZc3eY13a60/2r4mCiE24590G4Efwr9NYPBPusSu7+LLHNEWnucxPz/P448/TqFQIB6P9131vvGh9J04ueM2hfvAeEMAxJKGzI0u2RmP1NN9lOpwQ3dIGbO/ZynRFfyp5iDuOM4Y8C/A9Vpr9wLPeyPwCsITiE9orT+4wTG3AbcBaK1vrdXWqiksy8LzvNZ1pAsNQh9hMPrZrX3MfWWFoBKcd39sKEb2+WM7eq1u7WOr9WM/a7Uajz/+OKurqyilSCQS+L7f6Wa1VPk0LH4Dlr8F7moU2RWMXGOYfC6M3QyxxN7eo1qtcvXVV5PNZvfe4EgikQDY8lSk48HfcZw08Angj7TW/7iN5w8BXwZ+UWv90CaHmtOnTzduTE1NsbCwsIeWd79B6CMMRj+7tY/u/VWCpaBx5Q9gAkNsIoa9wyv/bu1jq/VTPyuVSmO6XiwWIxYLh3oymQyFQqHDrds7vwT5O2xWj9tUH29K608FZGdcDv2rBLVE6/pZq9U4cuQImUymZa95+PBh2Ebw7+igjOM4U8DfAu/QWn+p6f448BngfVrrf1r3tPrpZX5/WimEqItfbmMWq5jANFL+yoT3i/5VLBaZm5ujUqn03Xi+8aH4UFStf78Ffhg3Y0OGkZtcsjMuQ08LUAoSmQS13j/HAfa34O8S4NPAYaDqOM4LgXHgSuDtjuMAfFhr/bdRu64FDjQ9/y3AiwBDeFJwGiHEvoqlYlgzyUa1fyylpNq/TzVP16tWq31XuV99KhYuwnOHhZ+Pvr/KMHyVR/ZZLpmjHrE+Pqfd97T/PpK0f58ahH5KH/tHr/Vzo+l6W+mVtL9fUOTuDNfWrz6xltZPTPtkZzxGbnWxRy8cE1vdz4FN+wshhOgOQRA0dtfrp+l6xoPigxarsxbFb1sQRGn9lGHk5iitf2nQc9X6e9X7/7JCCCF2zfM8zp49y8rKCkBfTNczBqpPxsgdt8nfYeEXo7R+zJC+1iM745K+rr/T+luR4C+EEAPIdV3OnDnTV7vreXlF/g6L1Vmb2lNNaf2L/HARnls9rGx3DHUb1xDM+7i5Gu5odd9rZyT4CyHEAFk/Xa/XU/uBB8VvW+RmbYoPxtfS+sOGbJTWT17SXWl94xqCUx7GGIgZgqVg31fK7O1/dSGEENtS312vH6brGQPVJ6Jq/TttglIU2WOG9HVNaf0u7WKw4IcrZdK5lTK79KMRQgixV/02Xc/LKXInwmr92pmmtP4hn9EZl5FbPKyR7kjrb6oGSqlw4npExRSmvH9N6N1vgRBCiA0FQcDy8vI50/V6dUvdwIXifRarx21KD8XBhFfL8XTAyC3hVf7QxecvOd3VEmHqv5kJDLHU/o1NSPAXQog+4ft+Y7qeMaZn0/vGQOWxMK2fv8smKK+l9TNHXbLHPNLXeKje6xoAsak4QSEa84eOrJTZox+dEEKIOs/zmJubI5fLAfRs5b67qsgft1mdtXDPrvUheUlYrT9ys4eV6YG0/haUrYgdsVBzHiQUsdGYVPsLIYTYnmq1ytzcHIVCoWen6wUuFO4Nx/FL32lK648EZKO0fvJQj6X1t0HZith0HPtIAjuz/9thS/AXQogeUyqVmJ+fp1gs9mQRnzFQORUtwnOXTVCJqt7jhvT14fS89NU+qvfOZXpGb31jhBBigNUr9yuVSk8W8bnLitwJm9ysjbuwluJOXuozesxl5GaXeLqDDRwgEvyFEKKLGWNYWVlhYWEB13V7LugHtSitP2tTeqQprZ8NyN4aFu8lL+q/tH63k+AvhBBdKAgCFhcXWVpawvf9nkrvGwPl78bJzVoU7rYJqlFa34qq9Wdchq+UtH4n9cY3SQghBoTv+8zPz7O6ukoQBD0V9N0lRe64Te64jbu4ltYfeppPdsZl5EaX+HAHGygaeuMbJYQQfW79Rjv1P90uqMLiPTD/1RTlR9dCijUahJvpHHNJTPf+9Lx+I8FfCCE6qFwuMzc3R6lU6pmNdkwA5ZPxcBGeeyxMTQFWmNa/wVtL63f/ucvA6v5vmRBC9KH6RjvlcrlnUvu1hbW0vre8FtnTVxgyN1fJ3OgST3WwgWLbuv/bJoQQfWK7lfvGNQQLPtSARLgcrLI7syetX4HC3Ta54xblk01p/fG1tP7EkTSFgtuR9vUqzwuX9+3U0I4EfyGEaLOdVO439npX4c5vxjUEBS9cDnafTgBMAKVH4uSO2xTusTBuVK2fMIw80yN7zCX1DEnr70QQBPi+TzweZ2hoiMnJSUZHRzuW8ZHgL4QQbeL7Po899hiPPfbYtjfaaez1rqKAqxSGMBMQP9Te/7JrZ5vS+itrkT31jDDgjzzTIzbU1ib0DWMMnuehlCKZTDIyMsLY2BhDQ0ONf9tO2rfg7zjOEPAq4E3A7VrrjzqO8yzgfUAc+IbW+s2bPH/bxwohRCe5rtvYaGdkZGRnqd0a5wUHpVQ4BNAGfhnyd4UBv3JqbeK9NdFUrT8p1frb4fs+QRBg2zbDw8OMjo6SyWS6ctbGfrboIFAG/qbpvg8Dr9Bafx9wk+M4z9vk+Ts5Vggh9l25XObUqVM8/PDDFAoFLMva+WY7CRpbvdYZYyDRunaaAIoPxnnqY0OcvD3D/N8PUTkVRyUM2RmXS15b4ulvKTL14poE/k0YY3BdtxHwp6amuPLKK7nqqqu49NJLyWazXRn4AdT6L1m7OY5zO3AK+O/AF4AXA+8Frgc+rrX+wAbPGd/OsY7j3AbcBqC1vrVWWztVtiwLz/Na36EuMgh9hMHop/Sxt6ysrHD69GlKpRK2bZ9z5R6Px/F9f9uvZWoB5YdLoFSY8jcGjCF15TAqsbdAUn4Klr4BS98Cd2WtjSNXGyaeB2M3Q3wXG8zttI+9Kh6PU6lUABgaGiKdTjM1NUU6ne6KVD5AIpEA2LIxnR7zPwx8DHgL8K/YPBOx5bFa6w8TZggAzMLCQuOxqakpmm/3o0HoIwxGP6WP3W+jyn2A5osOgEwmQ6FQ2NlrHzYEC9451f7FWmlXqX+/FKX1Z20qj61lIezJgOyMS/ZWF3sivAgsu8AuivZ308deUS/UsyyL6enpxvh9PaNTqVQaJwTd4PDhw9s6rmPBX2u97DhOAXgl8BTwbuB3ABzHiQOfAd6ntf6nzY4VQoj9tB9r7itb7am4z/hQ/E64CE/xPgvjhxeCsaQhc5PL6IzH0BGfLrlY7Sr1Qr1YLEYymSSdTjM2NkYymez5E9Zm+1nwdwnwacIr+KrjOC8EXgN8ivBc8wta69mmdl0LHGh6iQsdK4QQbed5HvPz8+Ryua5dc7/6VCys1j9h4eej5KgyDF8VrrqXOeoRa2HtQL+oz7lPJBIMDw8zNjbG8PBw147Xt8K+j/nvI3P69OnGjX46Y7uQQegjDEY/pY/do1qtMjc3R6FQ2NV6++1OiftFyN0ZpvWrTzSl9Q8EjM64jNziYo+39//5Xkv7r59zPzIyQjab3XKr5F74zkZp/64f8xdCiK5ULBY5e/YsxWKx667yjQ/FB6JFeL5tQT2tP2QYuTncMnfoskDS+hFjTKMgsRvn3HdC93ybhRCiw4wx5HI5zp49S61WIx6Pb3k1uJ+qp2Osztrk77DwC01p/Ws8Ro+5pI96xLqnuR3VPOc+lUo15tzveOpln5LgL4QYeEEQsLy8zOLiIp7nddWVvpdX5O+0wrT+6bXAlTjoh9X6t3hYo307fLtt6wv1xsbGGBsbq099E+t0x7dbCCE6wPd9zp49y8rKyraX390PxoPCAxa5WYviAxYEUVo/Fab1R2dckpdKWr++TkQikWBoaGggCvVapfPfciGE2Gf15Xfz+TxAV6SCjYHqk2vV+kEpCmAxQ/rasFo/fb1HbID/126ec59MJju+OU4vk09MCDEwyuUyc3NzlEolYrFYVwR9L6fI3WGRO25Te6oprX9RU1o/O5hpfSnUax8J/kKIvlcv4qtUKl0xnh94ULw/DPjFB+Nraf3hgOzNHtlnuSQvHsy0vhTq7Y8tfwMcxxkjXFnvs1rr+zd4PEG44t7btNal1jdRCCF2rl7Et7S0RK1Ww7btjlbuGwPVx6Nq/TttgnIU2WOG9PUu2RmPzLUeasAuyeqFevF4XAr19tF2vmaTwC8CL3Yc57jW+i31B6JV+/5f4BEJ/EKIbrBREV8ng35tBZa+nCB33KI2t3b1mjwcpvVHbvawRgYrrT+IK+p1m+2eY35ba/3zjuMox3GSwAzwCuBG4P1a679vWwuFEGIbuqmIL3CheJ/F6nGb0kOACbfKi2cCRm7xyB5zGbo46Fj79tv6FfUmJye3taKeaJ9tJ5gcx/kscCVwGeHeUu8EXqe17v99HIUQXatbiviMgcpjMXKzNvm71tL6Km5IXxeuupe+xkcNwNB1PZWvlJJCvS61afB3HOcbwDyQAV6itS46jjMKXA/8KPB1x3HeqbX+7+1vqhBChIwx5PP5rijic1cU+RM2q7MW7tmmtP4lYVr/0PcnqdA9W762S/M0vPqceynU615b/ba8GPh1IA3MOo7zV1rrdzuOs6K1/m3Hcf4L8FeO4zxXa/2f2t5aIcRA22glvk6kjoMaFO4LV90rPRwHE17NxkcCsreGaf3koTCtb2WS0Dt73mxb89W9bduMj483tr4V3W/T4K+1zjmO85Na62dFV/xfdRznGPBnjuP8EPA7wF8C/8ZxnINa67l9aLMQYsD4vs/8/Dyrq6sdW4nPGKicirM6a1G42yaoNKX1j0Zp/av6O61fn4aXSCRIpVKMjY2RTqeZnp7u+t3uxLl28tvzg8AngacRbhf4fcAJrfWngE+1oW1CiAG31+10W8FdVuGqe8dt3IW19x+6zCd7zGXkZpf48L43a180r5c/NDTUuLqXQr3et9WY/6eAZziO8w/A0LqH3xAd4wCf0Vp/qD1NFEIMmmKxyPz8PKVSqSPj+UEVCvdarM7alB9tSutnA7K3hnPykwf7s1q/vl6+bduk0+nGNDwp1OsvW/1G/Tbwt9Hf/9B0vwHeDLjAKPA/AQn+QohdM8awsrLCwsICtVpt38fzTQDl78bDav17LEw1SutbhkyU1h++ykf12VR0WS9/MG015n+f4zgVrfX9juOsP829FfgF4I3AV9vVQCFEfwnKAf5Jl1xsBTeowpEYS4UlVlZWGkFoP4O+u9iU1l9qSus/LVqE5yaXeGrfmrMvXNeVaXgDbjundndHf/8WYZC3gOdqrf/KcZx7gF/XWr+sXQ0UQvSPoBzgzVYxCqqpKk8sPknxZIHY0yysof1L7wdVyN8drq1ffnTtPa3RgOwxl+wxl8R0/6y6t34a3sGDB2Ua3oDbzm/aJx3H+TXgp4EDWuuPAK8H0FrfBbzWcZxJrfViG9sphOgD/kmXol9msbqAW/YJfJ9YPI5aBg61971NAOVH4+SOR2n9WpTWtw2ZG8Itc4ev6I+0/vrd8LLZLOPj4zINTzRsVfCXBKaACeCthHP684QL/9T9C/B5x3F+Rmv9vba1VAjRs4wxrK6ucub0Gaq1CraySdpJan50dV1r33vXFhS5WZvcCRtveS2ypy73yB7zyNzoEl9fztyD1u+GV19kR9bLFxvZ6sr/k8DXCCv9H4z+fjfwlehxQ7j6X0ECvxBiPd/3WVxcZHl5OZwyllDYrn3O2LIxBpVo7VizX4HCXTa54xbl7zal9ceb0vpTvZ3Wb56GJ7vhiZ26YPB3HOcQ8GXgi8CvAu8BNPDTWutfbTruvwJ/ups3dxznZURDCFFbDmqtr7zAsRXgm9HND0QHh5S0AAAgAElEQVTrCwghutBGm+zYto2ZMgQFD0MYeI0xKAOxqb2PPZsASg+Haf3CvRbGjdL6CcPIM8NV91LP6O20fvPVfX03vHQ6LVf3YseUMRc++3UcxwIcwlX8PgK8DrhHa31D9HgCuA+4Xmvt7qUhjuO8Cjiktf69Czx+Smt9ZIvXuA24DUBrfWuttpZLtCyrMX+1Xw1CH2Ew+tmrfczn8zz55JMUCgUsy9owKJlaQO1MDTzAgsRFCVRi98GrcgYWvwFL3woX5KnLXGWYfC6M3UJH0/rxeLwx/r5Txhhc1yUej5NKpchms0xOTnbd2H2vfl93qhf6GWV+tkylXTD4O45zM/AJYBY4rLV+QXT/vcC/JywATAIlrfVv7KWxjuMMAd8Cnq+1zl3gmEeA08AC8Gat9cktXtacPn26cWNqaqrvl58chD7CYPSzl/pYn5+/uLhItVrFsqxz0/quIVjww3H9RHiVr2xFJpOhUNjdovd+GfJ32eRmbSrfW8sa2JNraX17ojvS+jvt50ZL6Hb7Xve99H3di17o5+HDh2Ebwf+CaX+t9Z3AFQCO48yue/hh4LPADPByx3GyFwra2/Q64O82ew2tdb0tLwXeAfz8Ht5PCLFHvu9z9uxZVldXLzg/37iG4JSHUaCUCm8XPGJHdj6lzwRQ+k64CE/hPgvjRWn9pGHkxnDVvdTTfXptqnrz2H0qlSKTyTA6OipL6Iq22mzMXwHXaa3vB5pPOY3W+gngCeCzjuO4wG8Ab99NAxzHyQKvJDyRqN/3MuDntNY/u8FTFLC8m/cSQuxdtVplfn6efD6PUmrTTXaCBb8R+CE6ASDKBIxv8/3OxMgdt8idsPFz0X9FyjB8ZTg9L3PUI9ZdWfAteZ6HMYZkMtnY/jadTssiO2LfbHb6PQl83HGcdwL/wXGcw2w8IedDhDMCdhX8gTcB/01r3ZwXuwi4tn7DcZxpwpkHVcJphq/b5XsJIXapUCgwPz9PuVze/nr7Nc4LaEqpLaf2+SXI3xmuuld5rCmtPxWl9W/tnrT+dtSv7uPxOENDQ0xOTpLNZuXqXnTMVgV/08BngEuBe4F7gKPA1dEhhjD9/mfAT0YZgW4hY/59qpf6WV/K1pRBpSB+uU0stfXYbbf0MQgClpaWWFpawvO8Ha/A5z/lYYrm/Kl9acXolWPnjIUbH4oPhWn94v0Wxg+fExsyjNzkkj3mMXSkd9L69av7iYkJfN9nfHy8bzfI6Zbva7v1Qj/3POYPoLWedxznxwiX9X2T1vrejY5zHOdOYJpwKEAIwblL2aqYIigbzGIVaya5rROATnJdt5HaN8ZsmtrfTGwq3pjap5TacGpf9alYuAjPHRZ+vimtf1VTWr8Hpq6vH7uvX90fOnSo6wOGGDxb/jZrrRccx/lF4NmEV/8bebXWurvnPwixz/yTbiPwQ/i3CQz+SZfY9d05SF0sFjl79iylUolYLLbnCnNlK2JHrEa1v0qo8ISgGmP+i3D2a8NUn2hK6x/wGZ3xGLnVxR7r/rR+c2V+87z7fry6F/1lW6fyWutZwil/DY7jJLTW9ZG7dziO806t9WqrGyhErzLltcBfp2IKU+5Qgy4gCAKWl5dZWlpqbKXbyg12lK2IH7LCtP4DFrnPWxS+bYGvgDixVJTWn3EZuizo6rR+89X90NAQExMTUpkvetKWv+GO4xwDhtfdbQO/DPyS4zijhPP+N1ycR4hBpVIQlM05JwAmMMRS3RHdPM9rTNWrp/bbEcQqT8bCzXTusPALa2n97FFD+qYK6aMesS6OnXJ1L/rRdk7v/5rwql8BP05YAJgBbMdxrgNuAE5orYtta6UQPSh+uY1ZrGIC00j5KxPe30mlUomzZ89SLBZbktrfiJdX5O+0yM3aVE+vpfUTB32yMy7ZWzzGLk5TKHTfaGFzZX4ikWB8fJyxsTG5uhd9Zbtp/18CcBznHq31LzmO8zDhev+/A1xOWO0vhGgSS8WwZpKNav9YSm272r/VgiBgZWWFpaWlxip8rUztAxgPCg9Y5GYtig9YEETV+inDyC0uozMuyUu6M60fBAFBEGBZFsPDw43K/G5eVU+IvdjOb/9GVTdGa/1Vx3E+QJgF+KvWNkuI/hBLxTpa3Oe6LgsLC6yurjaCWyuvYI2B6pNr1fpBKQqWMUP6urBaP32dR6y15xl71rzf/dDQEJlMRnbEEwNl019Jx3F+hnD9/o0euwm4GHix1jpoQ9uEELu0UdV+K69ivZwid0eY1q+daUrrH/IZnXEZucXDGumuav0gCBrLEKdSKUZHRxkZGZGrezGQNlve9wbCor7UBg9fBHwUeBnwsmjHPan0F6KD2l21H7hQ/HYY8IsPxRtp/Xg6YOSWcMvc5MXdldav78CWTCYZGRlhYmKi63bEE6ITNtvY517gJY7jPLDBw3ngeVrrkuM4Lwd+FPi7NrVRCLGJWq3WWJAHaGnVvjFQeTxM6+fvtAnKUWSPGdJHXUaPeaSv9VBdktavp/PrU/EOHDhANpslHo9v/WQhBshux/xLWutS9PPnCJf4leAvxD4xxpDP51lYWKBcLhOPx1sa4LxVRe6ETe64RW1u7XWTF/tkj0Vp/Ux3pPWbp+Kl02nGx8dJpVIyFU+ITWwn+CvHcV5LONVvMvp5NPobIAv8yLpFf4QQbeD7PgsLC6ysrFxwG93dClwo3mexetym9FAcTJTWz4Rp/dEZl+Thzpf3NC+0k0wmZSqeELuwneD/J8CB6OcPRz//adN9NnASuBK4v9UNFEKcOzd/q210d8IYqHwvSuvfZRNUoqvluCFzvUv2mEv6Gh/V4az5+mK98fFx0um0FOsJsUvb+d/jKq31G9bf6TjObwH3aK0/D7yt5S0TYsC1s4DPXa6n9W3cs2sBNHmpz+gxl5GbXeLplrzVrkmxnhDts53/SX4eOCf4RzMB3go8rx2NEmKQVSoVHn30UepbUreqgC+oQeHesFq/9EhTWn8kIHtrOCc/eVHn0vob7Yo3NjYmxXpCtMG2xvybbziOEwP+Anir1nqjmQBCiB0yxpyzAl+rKtSNgcp346zOWhTutgmq0Q6D8ahaf8Zl+KrOpfXr6Xzbtkmn07JuvhD7ZDfV/m8AVrXW/7UN7RFioNRqNc6ePUs+n2+swGdZ1p7Hst2lKK0/a+Murr3W0GXh2vojN7nE12/XtU/q6fx4PM7Y2BgTExOysp4Q+2w7wT/uOM641no5un018Ko2tkmIvmaMIZfLsbi4SKVSadkKfEE1TOuvztqUH1n71bZGA0ZuDefkJw7uf1q/PvdeKUUqlWrMvT948CALCwv73h4hxPaCvw1813GcdwNfAJ4CfsVxnPMO1FrLtr5CXIDruo3FeJqv8vfCBFA+GQ+3zL3bwtSitL5lyNwQjuMPX+mj9rkovjmdX98oR9L5QnSP7fzPUwVeAXyIcBe/TwKn2tgmIfpG81V+fTGeVlzl1xYVueNhtb631JTWP+KRnfEYudElvtHC3G3k+z7GGBKJhFTnC9HltjXmr7X+vOM4Rwl373sB8Pta6++0t2lC9K6NxvL3WrEfVCB/T1itXz7ZlNYfC8geC+fkJw7s36p7zen8ZDLJxMQEY2NjLd8qWAjRetv+LdVaF4CXRun/LzmO88JWVfs7jlMBvhnd/IDW+lMbHBMHPggcBeLAa6L9B4ToCkEQsLq6ytLSEpVKpSVX+SaA8qNxVmdtCvc2pfVtQ+aZ4WY6w1fsX1q/Ph0vHo+TSqUYGxuTnfGE6EE7PkXXWv+naLrfpx3HuTU6KdirM1rrH9jimJcBMa31CxzH+X7g/cCLmg9wHOc24LaonUxNTTUesyzrnNv9aBD6CN3Xz1KpxFNPPUUulyMIAlKpFKnU3nLu7kKM3NcyLH0Taktr4+SZKw0Tz4XxWyCestjFrzAAphZQO1PD1AwqoUhclEAlNg7gQRDgeV4jnT81NUUmk9nz+H23/Tu2yyD0cxD6CP3Vzx3P8wfQWv8Hx3EuA/4bcH7l3855juN8GVgA3qy1PrnBMc8DPuM4zkuA1xDOOljfrg8TLkEMYJoriaempvq+sngQ+gjd0U/f91leXmZlZYVqtYplWY1gWK1Wd/eaZcjfbZObtaicUtR/9ayJprT+ZJjWL/vALk+7jWsITnkYBUopTMFQWSwTO2KhbNXoXxAEJJNJ0uk0Bw4caIzfV6vVXfexWTf8O+6HQejnIPQReqOfhw8f3tZx2wn+f3qB+18J/LnjOJm9Xv1rra8AcBznpcA7CFcV3MjbgM8CLwEe3Mt7CrFTxhiKxSKLi4sUi0WAPY/lmwBKD8fJ1dP6Xhh8Y8korT/jknp6a9P6wYLfCPwQnQBgqM1XiR+0GBoaks1yhOhzWwZ/rfXtF7i/QDgLoJUUsAzgOM7LgJ/TWv9s9NgJwNZa/57jOM8FZHVBsS9c12VhYYF8Pt8Y795rUVttLsbqcYv8CRtvdS2yp64Id887+Jwhyl5lr02/wJtHAd8YPOMRUzGS8SST9hRTV0/JcrpCDICOl+U6jjNNOH2wCswDr4seugi4tunQjwMvchznS4AHvHo/2ykGy4WK9/YS9P0S5O8KV92rPLYWYO3JtbS+PRGm9eND7DqtvxljDF7MRQUxUlaK6cQ0GTuNMorYeEwCvxADQhmzf1OD9pmpb4wCvTFWs1eD0Edobz+b0/rGmL0vwuND6TthtX7x/nVp/RtdsjNemNZfV1mTyWQoFFoT/esV+pZlMTw8zGhqlOT9FsQUKqYwgUEZsGaSxFL7V7Uv39f+MQh9hN7oZzTmv2U1bsev/IXotOa0vuu6WJa15yvg6pkYuVmb3AkLPx8FVGUYvjIcx88c9Yi1cf2b9RvmTExMkEqlGuP8wbMC/JMupgyxlCJ+ub2vgV8I0VkS/MVACoKgUa3fnNbfS4GbX4T8nTarx22qjzel9acCsjMu2Vtd7PH2ZdqaA342m2ViYoKhoaENj42lYsSul9X3hBhUEvzFwGiu1i+VSkC4s9yeqvV9KD4Yrq1fuN8CP0rrDxlGbgrH8YeOBOel9VulHvATiQSjo6NMTk7KDnlCiC1J8Bd9r1KpsLCwQKFQIAgC4vH43tP6p2Ph2vonLPzCurT+s6O0fptmydUDfjKZbFzhS8AXQuyEBH/Rl1zXZXFxkXw+T61Ww7KsPS+16xcUuTvDtfWrTzal9cd8Rq6pkrmyip025yyW0yrrA/7k5KTMwRdC7JoEf9E3fN9nZWWlpeP4xofiAxarsxbFB5rS+ilD5qoqmWdUGTq4ltY3JlxEJ35o779azQF/dHSUiYkJCfhCiJaQ4C96WhAE5PN5lpaWKJfLwN5X3QOoPBlW6+fvPDetn7423Ewnfb2HOe2FK040UUpBbffvWw/48Xic8fFxCfhCiLaQ4C96jjGmEfBLpVJjPv5e5+R7eUX+DovVWZvaU2tp/cRFfrgIz60eVnatWt9PhOvkN29wY0y4Uc5ONBftjY2NMTExweHDh7t+PrEQondJ8Bc9wRhDqVRiaWmJxx9/nHw+35LCvcCD4rfDcfzig3EIorT+sCF7s0t2xiV5ycbV+rGpOEHBw2Aay+UqE96/5fs2TcvrZJV+UF6b769SyHx/IQaEBH/R1crlcmPFvfoqdSMjI3u6yjcGqk9Ei/DcaROUosgeM6Svd8ke80hf5xHb4i2UrYgdsQgW/HC9/IQiNhW/YLFffaU927YZGRlhcnLygvPw90NQDvBmq+EmPzFFUDaYxeq+r/QnhNh/EvxF16lUKiwtLVEoFBor7iml9jz27eUUuRMWueM2tTNNaf1DPqMzLiO3eFgjO1uER9lq0+K+esCPx+NkMhkmJydJpVK77kMr+SfdRuAHGkv9+iddWQBIiD4nwV90hVqt1pia57puSyr1AQIXiveH4/ilh+JgwkAXTweM3BIutTt0cdCKLpzD8zyUUgwPDzM5OUk6nT6nNqAbmPJa4K9TMYUpd6hBQoh9I8FftNx2x5HrAb9QKFCr1bYd8E0twH/KC6vqE5yXajcGKo+Fi/Dk77QJymtp/czRKK1/jYdq8bff88LS/+HhYQ4ePMjIyMie1hVoN5WCoGzOOQEwgSGW6q6TFCFE60nwFy211ThyrVZjaWmpsfjOTq/wjWsof6/UqLI3riEoeMSOWHilGPnjNrnjFrX5tbR+8mI/XFv/Zo94prVr6wdBQBAEDA0NMTk5ydjYWM9sixu/3MYsVsNd/Zp294tfLlMLheh3EvxFywTlgNpXypiVAGyFGo8Rs2NU3Spzd52lMlXdVcA/5z0WfOLKbqTQja8onLTJfy5J+XvWWlp/JCAbpfWTh1qb1m+emtfLc/FjqRjWTFJ29xNiAEnwFy1Rv+IPVgKUr6i6VVYKK5TGK3jKI16LY08m9h4ka4BSVJ6Kk3soSfGRBEEtKliLR9X6My7pq31UCy/A64V7lmWRyWSYmprqaKV+q8jufkIMJgn+oiX8ky5lv8KSWqTol/FwiREnXophZSxUcu/jyO6yYvWOIfL3JXFX1q5Ok9MeIzfUGH2BTzy957c5h+/7AF1duCeEEDslwV/smjGGSqXC4uIiudM53FoNa8hCVcHCRsUAn20vfLORoAqF+8JFeEqPNFXrDweMXF0jc1WF5HjQ0s106mn9+jj++Ph4VxfuCSHETknw3yZZCS1kjKFYLLK8vEypVGqkwlUCbNdGWQozZqAUYPxwHHmngdkYKJ+Mkztukb/bxlSjtL4VVutPPzeGSeVQHlG1/94Df/MCPKOjo0xNTfXkOL4QQmyHBP9tGPSV0IIgoFAosDS3RPGJAn7Nx0paxA+sbaBzzlK3cQWZGDHDjgK/u6TIHbfJHbdxF9c+16GnhdX6Ize5xFOQyWQoFFrz1a2n9dPpdCOtL4QQ/U6C/zYM4kpovu+Ty+VYXV2lXC7j13xij4d9t5WNKRmCU14juO90qdu6oAr5u8NV98qPrn0drdEg3EznmEtiuvXT83zfJ5VKceDAAUZHRyWtL4QYKB0P/o7jDAMfBi4GhoFf0lo/dIFjK8A3o5sf0Fp/aj/aOCgroXmex8rKCrlcjkqlglJqbfOcFYOJre1gp5TCYM7Zu36rpW7rTBCl9Wdt8vdYmFpTWv+GcHre8JV+WDPQIpLWF0KINR0P/lrrkuM4v6u1ftRxnNuA1wJvuMDhZ7TWP7B/rQv180po9UV3isUi1Wq1EfDP2zinxnlV7jvdu762sJbW95bXInvq6R7ZYx6ZG8O0fiv5vo8x5py0vlTrCyEGXceDP4DW+tHox8PAyU0O9RzH+TKwALxZa33OsdHJw23RazI1NdV4zLKsc27vhD/jkfvSCsRUI+VPYMjOjBEf7oqPENheH+tb4549e5ZcLofrusRiMVKp1KYbzlRHKvgF/7y96+OZOMnMhee7+2VYPgFL34DCI2vPtccNk8+FiefA0ME4EAe2N4RS3yRnsz56nkcymWRycpLp6emeWXWvbi/f114xCH2EwejnIPQR+qufypjWjqfuluM4PwX8KvAzWmtvi2NfCvys1vrnNznMnD59unFjamqKhYWFXbevF6r9L9THIAjOqdD3fb+xU952GTcc4zeKc/eu36CgzwRQeiRM6xfutTBulNZPhGn90WMuqSt2n9YPC/4K591fL96rL8LTLbvn7cZev6+9YBD6CIPRz0HoI/RGPw8fPgyw5X/uXXHZ6jjOvwVeSBjQvei+lwE/p7X+2Q2eooDlfWxiz62EVi/YW1lZoVKpEAQBlmURi8V2Vdy2nYK+2llFbtYmd8LGa1qEJ3V5OI4/8kyPWIsXxWueky/Fe0IIsT0dD/6O4xwjLPj7GvAFx3FqWusXARcB1zYdNw18EqgC88DrOtDcrlatVpmfnyefz1OtVoEwTdUo2tujjQr6/DLk77LJzdpUvrf2HtZEU7X+ZOuzS57nEYvFyGQyTE9Pk0gkWv4eQgjRr7om7d8GLU37dyNjDOVymeXlZYrFIolEgnK53PbxbRNA6TtxcsejtL63ltYfuTG8yk89vbXV+rB2lX/gwIFG1X6/XuX34/d1vUHoIwxGPwehj9Ab/eyptL/YviAIyOfzrKysUC6XCYKAeDyOUopEIkGttoPy+x2qzsXIzVrkTtj4uaa0/hUeozMumRs8Ym0YGalf5Y+MjHDgwAEOHz7c9b+AQgjRzST49wDXdVlZWaFQKFAuh4sL1LfFbfeVr19qSus/tpZRsCcDsjMu2Vtd7InWZ49kLF8IIdpHgn8Xqm+YU0/n12o1YrHYxvPv2/H+PhS/E1brF+9fS+vHkobMTS6jxzyGnu7Tjunyvh9OJ0wn0kzkx7CXE6gKMARsUbzfCzMyhBCiG0jw7xLr0/m+7zeu7vdrJbrqU7FwEZ4TFn4+CprKMHxVOI6fOeoRa0NdXX1e/tDQEFNTU2SHsgTHwyWViZlt7aUw6PsvCCHETkjw76B6Oj+fz1OpVID9S+fX+UXI3Rmm9atPNKX1DwSMzriM3OJij7enKNTzvPAqP51menqaoaFwHqB7f3XHeykM4v4LQgixWxL891FzdX6pVNr3dH6jHT4UH4wW4fm2BX6U1h8yjNzskp1xGbosaEtav36Vn0gkOHDgABMTE+fNTtjNXgqDsv+CEEK0ggT/NqsvtpPL5c6rzt/vjWWqp2Osztrk77DwC01p/Ws8sseitH6bmhQEAUEQMDw8zCWXXLLp1rm72Uth/XMCN8AsBZBQuPfL+L8QQjST4N8G1Wq1UaxXrVYxxmDb9r6m8+u8giJ/hxWm9U+vXWEnDvpkj3lkb3WxRtu31oPruliWxejoKAcOHNhWhiN+uY1ZrGIC00jfKxPev53nGN8QPOkDhviETbAUyPi/EEI0keDfAvW18+vFeq7rNsbu9zOdX2c8KDxgkZu1KD5gQRCl9VNhWn90xiV5aXvS+nBuAd/BgwfJZrM72kcglophzSQblfuxlNryyr35Od53PUhBfCI8WQgW/HBvgq8EJJ6fkhMAIcTAk+C/S81z76vVamOzHDyILShMLcBPBOetf98uxkDpMZj/cpLcCYugFAW4mCF9bVitn77eI9bGf/F6AV99yd1kcn8L7er7L5gyUDVhwD/tYYg2I1oJZwRIBkAIMegk+O9AqVRieXmZcrlMtVptFOvV0/nn7XznGoKCt+HOd63i5RS5Oyxyx21qTykgnIuXuMgPF+G5xcPK7iytb1zT2MCHBJuewNSv8m3bvmAB307tddpeffw/WPbXAr8xKFthFDIDQAgx8CT4b1MQBDz66KMkEokLFusFC34j8EMUdAgD6foNcfbUFg+K94cBv/hgvJHWj6cNIze5ZJ/lkrz4/LT+doL6dk9g6ivwDQ8Pc/HFF5PJZFrWv71O22uM/7vmnO2H1XhMZgAIIQQS/HekXqV/QTXOe1wpFQbbPTIGqo/HWD1uk7/DJihH7xMzpK93yc54HJwZolSpbvz87Qb1LU5g6ivwZTIZDh482JYZC3udtlcf/w++EmBWApStUOMxYnZsy1kDQggxCCT4t1KCxtVmnTEGldh9sPFWFbkTNrnjFrW5tXR68nCY1h+52cMaCdP6m43nbzsrscEJDIBbcbHUWmq/nbMWdjPVb71YKkbi+alzhg+2M2tACCEGgQT/FopNxQkKHoZz082xqZ2NgQcuFO+zWD1uU3ooTrjOLcQzASO3hDvoJQ8HO2vcdrMSTScwgQnwjc9QbIjD04cYu3JiZ++5S7uZ6reR3cwaEEKIQSDBv4WUrYgdsRrj6iqhtl3tbwxUHouRm7XJ32kTVKLnxA2Z68JV99LX+Kjd1tJtMysRm4pTy1VBwUgyy1RigoSysa7ZvwK5Vgbt+gwAIYQQayT4t5iy1Y6K+9wVRf6EzeqsjXt2LbglL4mq9W/yiGf2vgjPVlkJYwy+72MnbC665RDZ5QyqEjtvd7z92jlPgrYQQrSPBP8OCGpQuC9cda/0cFNafyQge2u41G7y0A7T+lu4UFbCxA1uySWZTzAdv4jMaJr4JTaxi88P6LJznhBC9AcJ/tsQlAPcR6t4T3qooc3nvV+IMVA5FWd11qJwl01QjQrv4ob0UZfsMZf01XtI629Dc1bC930CAtJ2mokzY9hxG2XUpkvhys55QgjRHyT4b6F+tRtgwANT3NnCPe5SvVrfxl1oSutf6jN6LNwyNz7czh6cy/M84vE4k5OTTE5OEjzoEcSDbQV02TlPCCH6gwT/LTSudnewcE9QhcK9FquzNuVHm9L62YDsrS7ZYx7Ji1qb1t/MZmvt+2Vv2wG9FVPwhBBCdJ4E/y3Ur3aNWSu622iKnAmg/N14WK1/j4Wpp/UtQ+ZoWK0/fJWP2seh8fo2uul0munpaVKp1HnH7CSgt2oKnhBCiM7qaPB3HCcOfBA4CsSB12it793tce1QD440xcLmKXLuoiJ3PErrL61F9qGnRYvw3OQSPz/mtlV9Fb7tbKO7k4Au8+aFEKI/dPrK/2VATGv9Asdxvh94P/Ci3R7nOM5twG0AWmumpqYaj1mWdc7t7fJnPHJfWiFQhkTZxlIWQdVQeSrN3P9QFL6zdlZgjxsmngOTz4Ghi2JAMvqzP3zfZ2RkhOnpaQ4cOLDtVfj8SY/ygyWCYkAsHSN1zTDx4U2+Gpe2qMG7tNt/y14ifewfg9DPQegj9Fc/VXM6e785jvMB4AtEV/PAdVrrp+32uHXM6dOnGzempqZYWFjYVTuDckDt4Qr3zH6P0ndSFB9OYNworW8bMkc9ss9yGb5if9P6sDaen0qluOqqq/A8b38b0AF7+bfsFdLH/jEI/RyEPkJv9PPw4cNwTq56Y52+8gd4G/BZ4CXAgy04ruXO3lvj3j/PUzk71rgv9XSP7IxH5kaX+NB+tiZUH8/PZDJMT08zNDTE2NhY138xhRBCdF6nB2tPAHdqrX8POAY8AOA4zsscx/nEVhXoZ/4AAAqgSURBVMftl+RYjMrZAGvMZ+KHqxx5S4FLX19m9Nn7H/h93ycIAkZHR7n66qu57LLLGBrqwNmHEEKIntXpK/+PAy9yHOdLgAe8Orr/IuDabRy3L0avtHnOu8eZMyexE535yDzPw7ZtpqenGR8fb+uuekIIIfpbR8f826xlY/4QptkffPDBTSvnW615PH96eppMJrPp8b0wHtUKg9BP6WP/GIR+DkIfoTf62Utj/mKdjcbzhRBCiFaR4N9FdjI/XwghhNgtiS5dwPM8LMviwIEDTExMyHi+EEKItpLg3yHN6+1fdNFFjIyMNNbbF0IIIdpJgv8+2856+0IIIUQ7SfDfJ77vA5DNZpmensa2ZTMcIYQQnSHBv808zyMejzM5Ocnk5CTxeLzTTRJCCDHgJPi3gTEG3/exbZtDhw4xNjYm4/lCCCG6hgT/FmpelOfiiy/eclEeIYQQohMk+LeALMojhBCil0jw3wPP8xqL8kxPT8uiPEIIIXqCRKtdqBfxHThwgMnJSVmURwghRE+R4L8D9Z31pIhPCCFEL5Pgv01KKa688kpZlEcIIUTPk3z1NimlJPALIYToCxL8hRBCiAEjwV8IIYQYMBL8hRBCiAEjwV8IIYQYMBL8hRBCiAEjwV8IIYQYMB2d5+84ziuBXwAmgA9prf9sk2M/CtwErAAPaa3/3b40UgghhOgznV7k54vAXwIZ4FHggsE/8kat9T9f6EHHcW4DbgP4P+3dfbBVVRnH8S8QCKSDZWASGpoNvdBEDSSRYSbJJKQ04i8ocgjfCKbQqdAcK1PyJXAKyklpKCqT/GGKNUwjk/g2Gs6EVkJkphVZmfiaiSAi/fGso4fjJQounO7ez2fmzL1n7333Wc9dF9Zea6/9LNsMHDhwu/2t76uoDjFCPeLMGKujDnHWIUaoTpx7Zdhf0pWSftnyervtB21vAwYA63dymkeBSyXdImlyRwfYXmh7uO3hQLfml6TVrduq9qpDjHWJM2OszqsOcdYhxi4W507tlZ7/fxqil3Qg8B1g2k7O8dly/P7A3ZJ+YvuZTi1oSimlVAPtvud/GNHwz7C9tmn7q4BlwEzba1p+bCuwGXhurxU0pZRSqpB23/O/oZThckkAF9heCfQB3gTs3zhQ0nzgbUTjP9P2lv/xsxZ2Son/v9UhRqhHnBljddQhzjrECBWKs9u2bdvaXYaUUkop7UX5nH9KKaVUM9n4p5RSSjWTjX9KKaVUM+2e8LdHSeoBLACGAj2AT9q+t72l6lySNgGrytv5wF+BeUS8v2g8ItmVSOoNnAZ8Bjjf9mJJ76Ilrq5cvzuI8XxgEvAw8JTtE7pyjACS+hKTpF4H9AVOBvpRrbrsKMbJVKwuASTdCPQCehMJ1fpQobqEDmM8kQrWZdV7/hOB7raPAs4BLmtzefaEh22/r7yuJ/4TmmL7PcAwSaPaXL5dcSDwLPDDpm0dxdWV67ejGAEuKXV5QnnflWPE9kbgS7aPJrJ5zqBidbmDGKFidQlge2yJ82fAKCpWl9BhjFDBuqz0bP/yeOAKypUZ8Bbbr29vqTqXpD8AfyMyIH4FuAIYC8wF3gossT2/fSXcdaUn/CfikdAVtMQFHEYXr99GjKXnfyYwhbgoWGp7QZX+hkus24DxVLAu4cUYnyCyrFWuLiWNA75MPIZ9FJGPpVJ12RLjaEBUsC6r3vMHOI9YEOh44Pk2l6XT2T7c9mjgKmAOMLB8v4CX9yq7sh3FVZn6tf31kpr6A8AUSUPKri4fo6QJwDuI9TsqWZdNMV5e1bq0vbzE9Wmio1G5umyNsap1WfXGfzVwj+0LgOHAujaXZ0/qBvwR+BdwKvAb4DjgznYWqjPYfoKO46pq/W4tX5+mAjFKOgU4CTjJ9gYqWJctMTY3BJWqyyZPAhupYF02eZLo7TdUqi6rPuzfE/gucDBxZTbd9v3tLVXnkTQA+DGR7vgRYCbR87gY2AKsKH+cXYqkQcRw4kAittuBxbTE1ZXrdwcxrgOOJYbGr7DtrhwjgKThwF3AHcALRFruS6hWXXYU481Ury4PJf4dbgL+CZwJDKFaddlRjCdTsbqEijf+KaWUUnq5qg/7p5RSSqlFNv4ppZRSzWTjn1JKKdVMNv4ppZRSzWTjn1JKKdVMNv4p1Zyk/pK6l+/fWb4OlNStvSVLKe0p2finVGMlV8RKIjUrwFJJvYCPA8sl9Wk6dqikZU3v+0laxU5I6inpdklry/PROzpugqQxuxxMSum/ls/5p1RTkl5DLF7yNdtXl20/An4K/J5IZ/oX4HwiX31vYuW6B8oplgFHAGc1nfaxkpGx8RndgW8TedKfJDocp9h+oemY6UTu9P2IpCm/BR4DehJJkAYTSaw2AqtsT++kX0FKtZU9/5RqqCyRfAdwCLBB0t2SVhON+yJgGrFc60HAxSW3+STgVtvDy/tXAvsSa0rMAa4jlilufEYf4AdEFrgpwBnAAODqsgxuw0FEzv/v2X617SOBQcSCKmOAPxOpY8cAs/bAryOl2smef0o1JOl0IiXtZcBU2w817buCSDk8m2i47wHuJzoLvYjUpyuB9wMjgGttj5X0c+As2/dKGkos9PIg8FHbz5Zz9yRGAt4NnG771qbPHE+s4gixaMqpRLrqRn71A4HP276q038hKdXMK9pdgJTS3md7IYCkXxMNfbMZxFB/Y98s4HctxzwHPEQs7HK4pB7A6tLwzwY+Rdw+OA5YJ2lf4sLh8fLzS4CFkq4jRiAGAx8s57uBWCgG4GPAa4mLiON3N+6UUsjGP6V6GwrcLKmjfb2ALwLLgQ+37LvP9kWS9iduHfSyfXbZdw3wTdsbGwdLOhUYavvMpm3nEbcNRhMLqCws5Zlj+/GmMo0BbtqtKFNK28nGP6V6G2J7MICkQ4A32r6pvF9MTLT7KjGxbzuSpgL9iBXdRgC3lV2jgLktFxR9gZ6SJrac5nPAeuAcYAPRuz9N0vXAM+WYfYhbDSmlTpIT/lJKDUOBCa0bbc8EvgF8AniUGA04xvYaYqj+ImJ4vnH8EtuDgDcA7y3fzwYWle+Ptj2ovJYQIwdvJiYbziPmGZzNSwYQs/9TSp0kG/+UUsMIoLekfZo3ShpGzNbfUjb1By6UdAAwFrgUGFZGDppdBMxtOVc/4DZJ210sAFtsj7Q9kphf0J2YH7CNGEl4Hnial0YDUkq7IRv/lOptkqQDJJ0LTCaepV8j6UrgWuKJgOnEPfn9gK3E43u3AN8nHgPcQvTYF0vqDSBpPPARYvLgi2w/RUwCXCBp2k7KNguYCtwHLAX+Yfv63Q04pZSNf0q1VLLuLSKer7+TuK8+0vYsYvh/JXE//g7gC8AFxMS/5SVBz4eAtbavA7C9lEj+c4mkI4DFwETbj5SP7E/M5Mf2PcQFwGWSji37e0haVTIGjiMyDt5ITAgcBxxJXKhcuId+JSnVSk74S6mGbG+R9C3gEdvrW/ZtJmbsXyPpUNsbgBNbTnGu7b+3bDuDmNjXG5DtVZLmEfMBNhEXDI3PuEvSMcC6smmc7cZjgEg6GHjA9q/Kps1E4z9416NOKTVkkp+UUkqpZnLYP6WUUqqZbPxTSimlmsnGP6WUUqqZbPxTSimlmsnGP6WUUqqZfwPA/CrYydrZSQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plot\n",
    "ax.fill_between(xs_pred, lower_95_percent, upper_95_percent, color='#888888', alpha=0.4)\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 预测区间估计 (y个别值的预测区间估计)\n",
    "\n",
    "statsmodels 提供了api直接计算出预测区间\n",
    "\n",
    "from statsmodels.sandbox.regression.predstd import wls_prediction_std\n",
    "\n",
    "$$\n",
    "\\hat{y_0} \\pm t_{\\alpha/2}s_e \\sqrt{1 + \\dfrac{1}{n} + \\dfrac{(x_0 - \\bar{x})^2}{\\sum_{i = 1}^{n}(x_i - \\bar{x})^2}}\n",
    "$$\n",
    "\n",
    "There is a 95 percent probability that the real value of y in the population for a given value of x lies within the prediction interval. \n",
    "\n",
    "预测区间{除了X的影响,还有随机因素}比置信区间{$E(Y) = E(\\beta_0) + E(\\beta_1X) + E(\\epsilon)$, 其中$E(\\epsilon)=0$} 要宽, 从公式看, 开根号下多了1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAEkCAYAAADZ6MOlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xm8ZFdZ6P3f2lW7hjN2n3QnpNOZe4gxXKOAIl4HBokCIggu4RoQESKCCM7CRQ2IgAwqXEXAl5eIaLwLGbzKRXgVRQVRRCFM6c4cQkKSzjmnz1DTHtb7x9q7eld1nbnmer6fz/n0OVV1qvauOl1PrbWe9TzKWosQQgghRps36AMQQgghxN5JQBdCCCHGgAR0IYQQYgxIQBdCCCHGgAR0IYQQYgxIQBdCCCHGgAR0IfpAKfXflFJqG7c7qJTyku+/Lfn30HZ+dxCUUucl/+aVUg9Pvr9gsEclxGSSgC5Ef/w28IzNbqCUOhf4BPDNyUXvV0oVgOcAH1FKldtuf5VS6sOZn+eVUp/p7mFverzHgE8qpc7HvZekx/I/lVLvaf8QopR6ilLq9zM/P1wpZbbxOHNKqRNKqY9vcbsXKaWu2vmZCDEelBSWEWJzSqlXAy/pcFURCIA4c9mN1tqXKqV+FvjZzOX7gCpQz1z2XGvtvyePcQD4KPB71to/Ty77C+CvgZPA9wNfA07hPhhcDZSAC4Dbkvv7MPAdwM9nHuMha+3STs95K0qpy4D/LzmHTyXB+9PAK5Pz/H7gbuDW5LLzgHlgLjkPgP8LnAu8OXPX37DWrmUepwR8CPccXAF8xlr7m23Hcj3wBOBA8hi3AV/GPRcAx4DbgRD4oLX2NXt/BoQYPvlBH4AQwy4JIO1B5AcAAxyz1p7q8GsHgD+w1v5Bp/tUSr0PmEq+/3bgT3FB/0Gl1H8CFqgA7wbegwvk5wNFa+0Lkt+7CnittfZpyc+vB2aA1yYPcyXwPuCNuzvzzpRSTwLeBuwHSpnjjYC/BG4EHgIuBJastU9Ofu8pwBOstS9Pfr4RmM0c76OAV+CeV5RS+4H3A0u4DynTwD8opWaBX7bWRsnvXQL8LjBlrX2fUmoe9yHge4EC8G/AY5NjXO7mcyHEMJGALsQOKaVmgHcCv7NBMN+pqwENvAX4qrX22zKP9Q6ggQt6x4HXK6XuB96Em+YuKKVuxk3VPw4XFP/SWnuNUurvcKP+bvtm4PHAx6y1fw+ka/0K+AhuZP4W4MnANUqpQ8DLgRyQSz4MfTg578cCb7LWPkcp9RXg75P7+m7ch5z/C7zUWhsDq0qpx+KC/KeVUtdZa79grX2eUupvgSNKqRclx/hTybFMA4vJ83AR8CPAv/TgORFi4CSgC7EDSdB6L24K+WNb3PyVmQDT7jDw/wBYa9+V3PcXaJ2SB3gxcD2QXve/cFPdP9t2uwZwD7CGC2w54HPW2i9ufVY7Y619U3K8/952uVVK/QhutPx53JLEi4HP4abes/LAHbhlg8PJaPzPrLUPKaXeBvwgLrg/CXiSUmofbmljJfn99wB/o5T6bdwHhSXch4xzgT9I7hvcB4YfwH3g+Y2uPAFCDCkJ6ELszOuBbwc+ADxFKfVC4ATw1swUcOp1W0y5t7sKN6Xc6VcKwG9Ya41Sahp4WNv1J6y1r0sC30VAwVr7q1udjFLqh4HXAX9rrf3FJLD+vLV2O8Hv25RSX9rgunngGdbaf08y4duP93PW2r9VSl2NW7ZYwiUOgpvO/wVrbZg5ztcCp6y12aS638HNUjwl+d3/g5t+f761tpF5Hp+J+0AkxFiTgC7ENiml3gj8BG5t9trk4tcD7wKepZR6gbX2psyv/KZS6uW4tfJZ4P7MdV8C/LaHOG6tvSR5rIuAo8mUNkqpG4AHktu9kTMJX9njex4ukAa4qfd/2sZpvQGXUPYapdQrcNPR280Un8kc71VA2Vr72eTnf8wc7ztw6+ntx3sNbvp+v1LqImvt3clVTwde1vbBZg6IlVK/1HY3zwJWgVcD/wG8FPgFpdT3ZW5TBGrbPCchRpZsWxNiC0qp2STj/OnA4621N6fXWWvvtNY+EbgB+GelVHZU/Gpr7RFclvUdgE5+/klcMP+JJIu7k6uAp3W6wlr7EtzU+0/ikuV+IzmuL+Gmql8H/Pg2T+/j1tqv46bGvwU3Xb2bLPDvSb46He/TcUH9Vbig+zvAo621D3U6Xmvtm6y1h5Pj+dbk+7cB1yfff5+19nDy9S+45Ytvwn04eCNuNuCVmUM4F5ekJ8RYk4AuxNbegZvNekQSNM9irX078Bjgpg7XreG2vX00yQj/WeBXrbXXWms3Gjk+CpdBXmy/IpmmvhY3Egc4CPyWUuoc4BpcwLw6GeVvylr7suTfurX2Wdbab7HWfn6r39vgeGeStfv2430q8OjMRVcCP6eUOorL3P9d3AzHXNuvvgtoWTZQSh0B/lUp9bjMOfwxsGitfbS19tG40XiAmyGYB47gPkAtcXaOghBjQwK6EFv7GWvtM621K5vdyFr7ZWttmlXuA7NKqR9SSv0hLjHrB3BT75cDVyXr3VnPUkqdo5R6JfBs3La1Lyml3onbDvZvye1ehMsan8VtFftT4B9xyXqvt9YGuL3dN2wyA9ANz1JKXaCUegtudHwR8NXk59fikvTALVP8z8zxvgW3L/x9wCustZXk+7erM1XyXgwcBX49+4DW2ltxhXY+lEzZb+bHcM/Dp4BP4vIEPru3UxZiiFlr5Uu+5GuHX7iAdf0m1/8HbkR4I/A/gP3J5Xng+bggU8dtP/Nx+83/E5dg92pgX3L7Ii4wfTK5/mBy+QdwI9AXJz+/B3hj2zH8MfD7PTj3/cBf4PIAbgJ+Abd+Dm4f/Atwme2fwCXnAXwmOd6nJj9/AnhJ5j5zwN/iPqg8DbgPuDxz/TuBF2R+fiou4/3hyc93Jo/xmeTyxyXf/xagcNPyHwN+etB/O/IlX736kkpxQuxCknUdWmuv3+D6hwO32I2n1EkKpOSttUtKqUcCD9gziWGdbn+ptfaODa4731p7X9tlHq7Yylqn39kLpdR3AbdZa7+xyW0us9bevsF1nY63hNuadhhYsNb+R1Ia9rtxe8kfn308pdR3Av9mrY2VUgvW2sXMdVcCDetG9OllHnCBtTatVCfEWJGALoQQQowBWUMXQgghxoAEdCGEEGIMSEAXQgghxoAEdCGEEGIMSEAXQgghxoAEdCGEEGIMSEAXQgghxoAEdCGEEGIMSEAXQgghxoAEdCGEEGIMSEAXQgghxoAEdCGEEGIMSEAXQgghxoAEdCGEEGIM5Ad9ADskvV6FEEJMIrXVDUYtoHPvvfc2vz9w4ACnTp0a4NH03iScI0zGeco5jo9JOM9JOEcYjfM8dOjQtm4nU+5CCCHEGJCALoQQQowBCehCCCHEGJCALoQQQowBCehCCCHEGJCALoQQQowBCehCCCHEGJCALoQQQowBCehCCCHEHllrCYKASqVCHMcDOYaRqxQnhBBCDANrLWEYEgRBSxCP4xjP6/94WQK6EEIIsU3pSDwMQ6IoQimF53nNAD6o0TlIQBdCCCE2lQ3iacD2PI9cLjfgI2slAV0IIYRos1EQH8RU+nZJQBdCCCEYzSCeJQFdCCHExIrjmEaj0QziSqnmuviokYAuhBBiosRxTBiGhGGI7/sEQTCyQTxLAroQQoixF8cxQRAQRRFRFOF5HkopcrkcSqlBH15XSEAXQggxltIgnk6nZ4P4OJKALoQQYmyka+JRFE1EEM+SgC6EEGKkTXIQz5KALoQQYuRkg7i1tpmdPmlBPEsCuhBCiJGQZqa3B/FxSWrbKwnoQgghhlYaxMMwBJAgvgkJ6EIIIYaGtZYoippbzKy1I1WtbZAkoAshhBioTkE8l8tJEN8hCehCCCH6Lu0lnq6Jw2jVTR9GEtCFEEL0xag3Pxl2EtCFEEL0zCQF8SiKWFpa4sCBA+Tz/Q+vEtCFEEJ0Vbb5ybgHcWsta2trLC0tUavVCIKAubk5pqam+n4sEtCFEELs2UbNT8YxiAPUajUWFxepVCpEUYTv++Tz+WY+wCBIQBdCCLErk1ZyNQgClpaWWFtbIwgC8vn8UM08SEAXQgixbZNWcjWKIk6fPs3q6iq1Wq25nc73/UEf2lkkoAshhNjUpJVcTdfFl5eXqVarAOTz+aEM4ll9C+ha6xLwQuAXgeuNMTdora8HngV8AzhtjPnhfh2PEEKIztLM9Gq12lKtbZyDOGy8Lj4q+nmk5wFV4M/aLn+DMeaGjX5Ja30dcB2AMYYDBw40r8vn8y0/j6NJOEeYjPOUcxwf43ieaaGXdDq90Wiwf//+QR9Wz0VRRKVSYXV1lUajge/7zM7O7vr+PM9jYWGBhYWFLh7l9ihrbV8fMBmV35mM0F8OXIsL9O83xrxti1+39957b/OHAwcOcOrUqZ4d6zCYhHOEyThPOcfxMS7nmQbxIAhatpcBLCwssLi4OMjD65nsurjnedTr9a4lttXrdS6//HL27dvXlfsDOHToEMCWUyMDTc0zxvy+MeaRwPcD12qtjw/yeIQQYtxZa2k0GlQqFdbX16nX68D47hNPxXHM6uoqd999N7fddhunTp0iiiIKhcLYnPewLA6kG/dWB3oUQggxhiap0Eu7arXK0tLSyK6L70Q/k+IOAx8GDgF1rfXjgK8CTwQs8GZjzL2b3IUQQoht2miP+CQE8WHfL94rfQvoxph7gEd2uOr1/ToGIYQYZ5NW6CUriiKWl5dZW1sb+v3ivTKe8w5CCDEhNtojPglBPI5j1tfXm3XU0/OepCCeJQFdCCFGiLWWKIqaddMnZY94VqVSYWlpiWq1ShzH5PP5sV0X3wl5BoQQYsil28vSkThMTlJbalLXxXdCAroQQgyhSeojvhFZF98ZCehCCDEk2reXpdPokxTE03XxxcVF6vX6xK+L74QEdCGEGKBJ3l6Wsta27Be31sq6+C7IsyWEEH02yZnpWUEQsLi4yPr6OmEYksvlJu456CYJ6EII0WOSmX5GGIbNdfFGo9HMC5DR+N7JMyiEED2wUeOTSZpKT6V11JeXl1vWxSWId5c8m0II0SWTXDO9nbW2uV+8VqtNxH5xG0Pj6z71BQvda7a2beP7zAohRB/EcdycSo+iaCKT2rLq9TqLi4tUKhXCMBz7/eJRVVG7xadywqd6skC85jH7YyHnXdr/Y5GALoQQOxRFEbVabSJrpncShmGz6Euj0SCfz6OUGsutZtZC474c1ZM+1RMF6nflwZ7Jg8jNh6gBnbcEdCGE2EKa1JZOp/u+TxRFEx3E4zjm9OnTrKysjP1+8bimqN7qU01G4dFKZrbBs5QuDSgfb1A+HhDPV7jgyOUDOU4J6EII4mpMdHuArYIqQ+4yH688nlOk27VZudVcLjdx2engnpNsMxRgLJPbrIXggZwL4CcK1O7MQ5wZhc/FLoAfCygfCfBKtnldvT6II3bG61UQQuxYXI0JP1vHKlCeIq5a7EN18o8qTlxQl6S2zmq1GktLS6yvrxPH8VjuF48bULvNBfDKCZ9oOXN+ylK8xI3Cp44H+A+LGMbPcxLQhZhw0e1BM5iD+9fGluj2AO+biwM+ut5Lk9rSID7pSW2pIAh44IEHxroZSnDK46HP5lj+4iy1232IzkRpbzqmfCxg6ooGpaMBubLd5J6GgwR0ISacrZ4J5inlKWx1QAfUB1KprbMoijh9+jSrq6t4nke9Xh+rZihxALXb/eZUerjYOgovXBgwddyNxAuHItSIfXaRgC7EhFNliKu2Jajb2OKVh3BOcZekUtvG0mYoy8vLVKvuU1w+n6dUKhEEwYCPbu+CRe/MWvjtPjbIjMLLMTPfZCkcqVA+GpCbGf5R+GYkoAsx4XKX+diH6tjYNqfblXWXjzJpP7q5tOhLtVodq6IvNoTaneko3Cd4sPWcCofCZkZ68cKQmdlp1tcbAzra7hr9V08IsSde2SP/qGIzy90rq5HNcpekts1lm6GM07p4uOxRPelTOVGgdquPbZwZhatiTPloQPl4QPlYg/zcaI/CNyMBXQiBV/ZGNgFO2o9uLtsMpV6vk8vlRn5d3EZQvytP5UTBjcLvbw1l/sNCyscDpo43KF4UoiYkNUICuhBipLSvhwOS1NZmo2YooxzEwxVF9aQL4NVbfGz9zAc2VbCUj5wp7pKfjwdyjDaw2IdigmqDYH+97zNdEtCFEENP1sO31t4MxVo70kVfbAz1r+WbCW2Ne9tG4ee6UXj5eIPSxSFqwKdpA0t8T4iNLeQt8WLc93oOo/lKCyHGnqyHb0+noi+j+hxFa4rqLS6AV0/6xNXMKNy3lC47Mwr3FwYzCt9IvBi5eg5qcPUcJKALIYZGdipd1sM3FoYhi4uLrK2tEYZhM4iP2vPk2o3mmtXZGl9vbXSSX4goX+GqsxUvDfCGecUgAMVg6zlIQBdCDIysh29fFEWsrKywsrJCrVZrBvFRm1Jvthu9OWk3up75EJKzlC5tMHWFG4n7B4ZrFL4pH2zYmkHf73oOo/WXIIQYebIevn1p0ZdsM5R8Pj9SyW3NdqNJRnr97rZ2o/uiZnW20uUBXmGAB7sH3kLOraFbF9QHUc9BAroQoufSeunpVDpIEN+ItZZqtTrSRV9a2o2eKBCtbtxu1D93OBud7JTyFd7hPOrBEAoKb78nWe5CiPGQrZcu6+Fbq9frzeS2KIpGal3cWgjuz1E9uUW70eMB5ctb242OE+Ur1Dke/uUF/H39r+sgAV0I0RXpVHq1Wj2rXrqsh3cWhiFLS0usra3RaDRGal08rrt2o2lxl+h05jX2RqPd6LgZ/r8aIcTQat9aViwWm4FcdJZNbhuloi/WQnjKawbw2h1t7UZnYsrHXAAflXaj40YCuhBiRzbbWiaBvDNrLWtra2d1NBv2kfhW7UaLFyY10q9oUDh/9NqNjpvh/msSQgycbC3bnfbktiiK8H1/6IN4sOhRvbnAqdt81k8uYMPMKHwqpnwsSWg7GpCbllH4MBnuvywhxEBIlbbdazQazeS2bEczFSmi+0MIAN9tc1L+4BeWm+1Gb/apnCwQPtj6Ia1wQUj5WIPyFQHFw6GMwoeYBHQhBCBdy/Yi7Wi2urraktyWros363wrV03Mhu5n73B+IEE9XPaonkjajd7W1m605NqN7v9vCnXxGvlZGYVvl7W2uQ99ECSgCzGh0qn0dGuZtVam0ndgJ8ltzTrfSWlQhcIqS7wYkTuv92/Du2k3Oj09zfq6BPPNZD/8+r5PoVBgYWGBqampgRyPBHQhJshGVdrSQC42l63c9uCDDwLbTG7rVOcb5abfe2TLdqNHk33hxwbXbnSUWGsJwxClFL7v4/s+5XKZmZkZCoVC8/9PGuAHoW8BXWtdAl4I/CJwvTHmBq31twNvBnLAvxpjfqlfxyPEpGifSpeM9J3plNy2b9++nSW3JXW+s0HdYrs63W6jpN3oyU3ajR5zGenD0G50mKWzV2kL2kKhQKlUYnp6mlKpNLQzWP18Sc8DqsCfZS57F/BUY8zdWuu/01o/xhjz6ewvaa2vA64DMMZw4MCB5nX5fL7l53E0CecIk3Ge/TrHdCSRjsSttZTL5b6MwPP5PAsLCz1/nH6o1Wo8+OCDrK6uEgQBvu8zOzsLuFmN6enpbd+XPRxTu6MKyUyItRaspXS4jCrs/oNVuAprX/VY/UqO9Zs9okpmLbxgmT4aM3tlzMyVMYVzLKCAYvK1uZ2e46jyPI9SqUQYhs06+aVSibm5OWZmZnZcHyCKIqampgZSV0D1ewFfa309cCfwV8DHgWuANwHfDNxojHnrJr9u77333uYPBw4c4NSpUz071mEwCecIk3GevTzHjbLS+21hYYHFxcW+P263pG1J19fXW5Lb2rn15fUd3bcN3Jr5XrLcXbvRPJVkX3jj67nWdqPnRM3qbHttN7qbcxwFcRwTRVFz6nxhYYEwDM+aOt/L/ZdKpa5uTzx06BDAlgc26EmXQ8D7gFcA3wPIHKAQ2yRZ6d0RRRGnT59mdXW1pS1pt0dYyle7SoCL1hXVW1wAr97it7YbzVvKl6Vr4SPWbrQP0tkqoGX0PTMzQ7FYJJfLjfyH0KyBBXRjzJLWeg14AXAf8AbgNwZ1PEIMOynw0j1xHLO6usry8jL1eh0YnrakW7Ubze+PzqyFXza67UZ7IfvhtlAoUCgUmJmZoVwuD31Bn27oZ1LcYeDDuFF5XWv9OOBngA/hJqE+boz5bL+OR4hRkJ1KT6cJJaFtd6y1VCqVZnKbtXZoyq/GtXQU7lM92dZuNGcpXZKMwo8H+Ael0QmcPXXu+z5TU1PNde9J3LXRt79kY8w9wCM7XPUd/ToGIUbBRm1HZRS+O7VajcXFRSqVCnEck8vlBv5cNtuNJsVd6ne1txuNkgDeoHwkwOt/J86h0j51XigUKBaLzMzMUCqV5ANuYvAfTYWYcOmbVbbAiwTxvQmCoJncli2/Osg3/i3bjV4aNLuVTXq70UmfOt8teWaEGIA4jls6loHUSt+rbG/xIAiaz+eg1sWtheBBr1ncZcN2o1cElI5MbrtRmTrvni0DutZ6Hy5x7aPGmC93uL6AS2h7lTGm0v1DFGL0pRXaarXaWWVWJYjvXjZDPVt+dVCjuLiRtht1QTxcams3elHQnEqfxHajMnXeW9v5qz8HeA5wjdb6P4wxr0ivSBLd/gS4VYK5EK3a94YXCoVmpTYZdexeNkO90WgMPLkteMhrBvDa7f7Z7UaPuoz0SWw3KlPn/bXdZ/Qrxphna62V1roIPAq4FvgW4C3GmL/s2REKMUI2SmjzPE/Ww/fAWsva2hqnT5+mVqsRxzH5fH4gz2kcuOpsSzdNUTlRIDzVod1oUtylMEHtRmXqfPC2/RFJa/1R4ChwEdAAXge8xBgT9ejYhBh6ktDWO51qqA8quS1c8qgkNdJrt/rYIDMKL8WUjgVMHXNJbbkJaDearXWeLdgy7LXOx92mAV1r/a/AA8AMrub6utZ6Hlem9QeBT2utX2eM+aveH6oQwyGt0BbHsSS09UCtVmNpaYlKpUIURc3Kbf18fm0ItbvyzYS29najpQtiikdrlI8HFC907UbHmUydj4atXolrgJcC08BntdbvNca8QWu9bIz5da317wHv1Vp/pzHm13p+tEIMQLZCWxzHktDWA0EQNDPUwzAkl8uhlOprsAhXVHMtvHprW7vRYkz5SNBsNzp/qMz6erVvx9ZP2eY+6fR5pzahYvhs+r/FGLOitf5hY8y3JyPzf9FaPxJ4u9b68bhSre8GfkhrfZ4x5v4+HLMQPbdRhTZJaOuedJvZ+vo69Xq9ORLvVxBvthtNG53c16Hd6BUuI7100Xi2G00/rALkcjl8329mnZ9//vmsrKwM+AjFTuzkT/SxwAeAi3FdX74L+Jwx5kO48q1CjKz0jU0qtPVWeyOUNFmwX3vFozVF9aQr7lI76RPXMqNw31I64gL41LGA/P7xa3SSJq6lU+e+7zM9Pd2x3adMpY+erdbQPwRcrrX+IFBqu/plyW008DfGmHf25hCF6I32bWWS0LZ7cS0mvitk1VshjBt4F+fxSi5YRlHEysoKKysrLXvF+xHEbQyNe5J2oycLNO5pfcvLH4iYOu7qpBcv2Vu70WGTHX237/kuFouyXNRF2fePjVru9sNWH8F+Hfjz5N8PZi63wC/hmqrMA38NSEAXQ22zUbi8ue1eXIuJPt/AKrDTMfF6TLRYo3Yk4HTVbTMD+rZXfLN2oypvKV12priLf874jMLbE9eKxSLT09OSuNYD2QCePt/pDoxB2moN/Uta65ox5sta6/a//EcAPw68HPiXXh2gEHsho/Dei+8KXXdPBSuNVR6sPkA1rMFJKJxf7HkwsXFbu9GvdWg3mgTwcWk3Knu++6s9gKe91AcdwNtt53/aF5J/fxkXuPPAdxpj3qu1vgl4qTHmmb06QCF2qr24i7Qc7R1rLeurFZary9TiGrlcznU083LQwwoVUVVRu9VvdiuL1zZqN9rAPxiPdKMT2fPdf6MSwNttJ6B/QGv9fODpwEFjzHuAnwUwxnweeLHW+hxjzEM9PE4hNpTWSY+iSIq79IG1tmWveFhroKIcnvLwcz6NqIHFovzuRVFrIfhGrpnQdla70Xk3Cp863qB0+Wi3G80mrvm+T6FQYHp6munpaZk675H2NXDf98nn8yM307FVUlwROAAsAK/E7TlfxRWbSf078DGt9TOMMXf17EiFyEj3hacBHKS4S6+1BPEwbK4Z+geLxPeEWNzrYLEoC97C3j5MxXWo3nqm0Um0cna70TShzT9vNNuNZpuVpMG7WCwyOztLoVCQv+ceSQtCpbN3xWJxJAN4u60+7n0A+BQuw/3m5N83AP+cXG9xVeTWJJiLXmofhQPNPeGj/p9wmNXr9WYQT/uKK6XIkyd+IHJpsT6oh+VgJQalUDmFt5Db8QjdtRvNJfvCfWp3trYbzc26dqPl46PbblQqrg1GNoC3j8CnpqaoVMajt9iGf0Fa6/OBfwI+Afwk8EbAAE83xvxk5nZ/BPxhj49TTCAZhQ9Go9FoFnwJgqC5dphuM7OBdSNyBQqFDS1UIrzDeUr7SkTr218837Ld6MWuMlv5igaFh41Wu9Hs6Dvd810ul5mdnZXEtR7bLICPsw0DujHmPq317wIaeD7wHuB64EfS2yS90B8P/FxvD1NMAhmFD05aenV9fZ1Go3FWEM+KF6NmMIckqCtLvBjBvm081inXbrRywqd2hw/ZdqPTsQvgx0av3WgYhlhrmwEkTVwrl8uSy9Fj7QG8UCg0ywdPks1G6N8KvB/4LPDvxpgXJ5ejtX4sLkmuCHzEGBP042DF+ElH4aurq6yvrwPDPwpPi6jYGqgSLUVURslOgnjrL54J5imFctPvHcQB1O/wXXGXEwXCh9rajR4OKR9rMHVFQOGC0Wg3GsdxS63zQqHA/v37mZmZmYiR4KCNSxJbt202Qv8v4AiA1vqzbVffAnwU1xf9WVrrOWOMFP0VW9osI32Yg3gqW0RFeYq4ZrFLDbi6MBJBPQxDFhcXqVQqLfXTd1S1zQcb2pag3p6P3V6HAAAgAElEQVTVvmm70XJM+ajbUlY+GoxEu9F09J1uG5ufn2dubk62jfWJBPDt2WyEroArjTFfBrLvVNYYcw9wD/BRrXUA/Dzw6p4eqRhZ47QvPC2iorxkutlT2NgS3xXiHR/OiiUbNUHZbelVbyGXrKG7oG6xEEJjucQ3Ppxn5UvzBA+0vrUUzg9dAB+BdqPZoi3p2nc6+k6fs4WFBRYXFwd8pOOr0z5wCeBb2yyt8hzgRq3164Bf0VofAhodbvdOXCa8BHQBjHd1Nls7E8xTylPY2oAOaAPtQbybTVCUr/AO52ncbandVqR2d4HaPYWz240eTRLajgXk54e3xGqauJZmnpdKJWZmZiiVSiP5oXMUjWohl2Gz2ZT7Ka31E4C/AS4EvgjcBNyrtb4juZkFng1UtNaHk5G7mDCTVCNdlSCu2ZagbmOLVxr8yKHbI/F2NoL63flmcZegvd3oeSFzV0H+0nVKFw9nu9HtjL5F7w1rLfRRt1Ut9we01k/ClXz9RWPMFzvdTmv9X8C5uGl4MQHGeRS+Ge/iPHapgY1tc7pdWXf5IKSJbXtaE99EtKqonkwy0m/Zut3o9PQ06+thVx67G9rXvmX0PRgSwPtjy3ehZKT+HOA7cKP0Tl5kjBme/8Wi6yZpFL4Zr+TB1YVmlrtXUn3Pct91dvo2tLQbPVGg8fW2dqMHI6bS4i6XBkM1Ck/3fXfKPJfRd39Za5vvERLA+2db/x2NMZ/FbV9r0loXjDHpmvpva61fZ4w53e0DFIOTbs1JM9LT4D3uo/CteCWv7wlwjUaD5eXlngTxlnajJ33iymi0G00/WKa5AeVyWRqWDIi1tln8KZfLNVvlSgDvry0Dutb6kcBU28U+8BPAc7XW87jCMq/p/uGJfkpHONlp9LSoi7xB9l+9Xm8G8U4V23ar2W705qTd6D0d2o1ekYzCLwvwhmBwm84QAdIudAhkA3j6NykBfPC2M0L/U9zoXAFPxiXJzQC+1vpK4OHA54wx6z07StEzaWGXOI6b1ZbSUbi8SfaXtbZZO71arXY1iG+r3egVrtlJ/sDg241mO46lNc9nZ2el6tqASAAfDdudcn8ugNb6JmPMc7XWt+Dqu/8GcBnw9t4douimNJmtU3lV+c/Zf9ZaKpUKp0+fplqtNruYpevAu79f1240XQuv3z287Ubb+30XCgXK5TIzMzMUi0X5YDkA6Wxdug6eJhXKe8Rw205A71TGyRpj/kVr/VbcaP293T0s0S2dktlGubDLOIjjmEqlwvLyMvV6vWtBfKt2o6XLgma3skG2G80mr6Wjb+n3PVhpElv63pDL5ZidnW2WtxWjYat+6M/A1WvvdN3VwAXANcaY4cmUERO7pWyYRVHE0tISX/va16jVasRxvOcg3tputEDtzvzZ7UaPN1xxl6MBXmkwJVY3Sl4rl8vyoXJA2qfQc7lcs5hLGsAnsbnJqNus9OvDcYlv5Q5XPwy4AXgm8Eyt9Z2S4T44myWzyRvm4IRhyOrqKisrK9Trdaampprr4rv9YBU3oHab3+xWFi13aDeaZKQXzu//KDz9W7TWNreOScvQ4dDekSyfz0vQHjObVYr7IvBUrfVXO1y9CjzGGFPRWj8L+EHgL3p0jKIDSWYbTuke8Wq1Sr1eB2hu4fF9n0ajU/XkLe5zy3ajbhq9fDQgN9XfUXj79Hk6dV4ul2X6fMCkpejk2e0aesUYU0m+/1tc+VcJ6D202TS6jMIHJ81MX15eplKptGSm7zagxQHUbvepnuzQblRZCocDptJR+AVRX9uNtmefp5XX0ulzaVoyONLQRGznHUdprV+M27Z2TvL9fPIvwBzwA22FZsQepclslUqF9fX1ia3MNozSpLbTp09Tq9W6ktQWLHlUb07Wwm/v0G70WJLQdiwgN9O/UXi2cUmxWGxOnxcKBQkUAybV2ES77QT0PwAOJt+/K/n+DzOX+cDtwFHgy7s5CK11DfhM8uNbjTEf2s39jLo4jmk0Gi3T6DMzM5LMNgSiKGJlZYXV1VUajQZRFO0piNsQanfmqZ50GelntRs9lGk3erg/7UbbS6dK45LhInvBxVa2E9CPGWNe1n6h1vqXgZuMMR8DXrXH4/iGMeb79ngfI0f2hA+3tNxq+3r4brf8BUuw+vmi21Z2q49tnBnhNtuNHncj8fzczkfhNrDEixEEgO/6lit/41F0p+1jafEWWf8evPatZBLAxVa287/22UBLQE8y4F8JPKZLxxFqrf8JOAX8kjHm9sxjXQdcB2CM4cCBA81fyufzLT8Pu/QNNK2PHsfxll2f8vk8CwsLfTzKwRiG87TWsra2xuLiIpVKpVkzvVQqUSqVdn5/EVTuVKx9OcfqVzzq97a+zsXzY2aujJm9MmLqMpuMwvNss95T62M1Ymr3VMkp33WBCy3cbyldWkYV3OOmHyCz0+f79u1jamqqazNAw/A69kMvzjMN4EBzJ0ShUGgutfXbqL2/7tY4nadKp3A2orV+0BhzMPOzB3wa+BNjzB9182C01j8C/Kgx5tkb3MTee++9zR8OHDjAqVOnunkIXZedRk8bnOwkE31SkowGdZ7p1rLV1VXq9XrL/vBd3V/SbrR6wqd6i4/NtBv1CpbiEVedrXwsIL+ve+UbovtDbNWiOHPcsY2JSzH5c1u3jxWLxZ6N8uTvdWeGeSvZKLy/dsMonOehQ4cAtvyj2E2W+8uA090O5gkFLPXgfvtms2l0WQcfvDQrfWVlpTkKB5pvpDt9jWwM9a/lm8VdGvdu0G70ioCFKwtUGz1qeRC4/6ihDcipHAWvQClXZKY0x8yRmaEIEKI1E122kolu205Az2mt9xtj0kB7HHhhtw5Aa30u8AGgDjwAvKRb990PUlp1+KW7BdKs9CiK9rS1LFprazdabWs3enngstLb2o16fgG6uA8k/XvL5/P4RZ+iLTJbnKXguQIuNrZ4M4OZrhWOZKKLftrOu5kP3KG1fgPwceA+4Hla67NuaIzZcQtVY8wDwHfv9PcGabOiLjIKHw6dEtrSIL7TN1MbQ+PrOaonXXGXRnu70YWomZHey3aj2RKqxWKR+fn5ZgU2e4El+nwDq2gGc2XBu1iS2/pJMtHFIG3nf3sduBZ4J6672geAO3t4TEOnvagLSDb6sElH4SsrK8294XsahVcVtVt8Kjf7VG/p0G70soYr7nKsd+1GNwrghULhrNuqkoKrC8R3hdgaeCWFd3EeryR/n73UnokuXcnEIG1rDd0Y8zGt9VW4rmrfC/yWMeZkbw9tcNpro0tRl+GTXQtvH4XvZm+4tdC4L+eqs321QP1rraPw3L6oWZ2tdHmAd3ZM3bM0cTJtIbpZAO/EK3l4x3twYKJFdmauU1MTIQZl20MXY8wa8CPJ1PsntdaPM8Z0qvM+ctrXwbNTZjKNPjzCMGR9fZ3V1dWurIXHNZW0G/WpniwQrWQ+rHmW4vkNShfVKB9uUHg4eIXuvmGnATxNjpqammJubk6amAyZTolss7Ozu6rLL0Qv7fhd0Bjza8nWtQ9rrR+RBPqRValUzqqNLm+mwyGOY6rVarNbWTYjfVdr4RaCB9rajcat7UZLF9UpHapROhzgFdwHO4vFLik4b2/r0XEcEwRBMzlqZmaGubk5isWi/M0NkXQdPJ2R65TIJq+XGEbbquXefoEx5le01hcB/y9wdnbciMhmoIrh0Gg0OH36dHMaPbsvfFej8K3ajV4SNHuGF86PiL8eQth6Hwrlqq/tULYSW7oGvm/fPkqlkgSEISIV2cS42M475B9ucPkLgHdorWdGfZQuBiedRj99+jSnTp1qSWbbbd/wlnajt/sQZRqdzMSUjzaYuiKgdDQgV24rs+CDDVsLtFjspiVUm7dLAjhAsVikWCwyOzvL1NRUXzuRxbW4mRynSkhyXJv2TPS0GqAsrYlRt2VAN8Zcv8Hla7jsdyG2Le1Ulq6DB4Eb+s7Pz++60Umz3WiyFn5Wu9ELg2ZG+lbtRr2FHPE9IVa5oG5Jtn8tdH6zTzPRfd+nVCoxNzfX1VKqOxXX4jPb1zxFXLPYpQZcXZjooN6+Dj5MFdmE6BbZpCp6ylpLrVZrbier1+tYa/c0jQ4QLHpn1sI3ajd6vEH56M7ajSpf4R3ON5ucKF+1NDnJJrIVi0X27dvX3As+DOK7wmYwB/evjS3xXeFEZcC3F3SR3uBiEkhAF9uy3Wnc9u1kaR37rbLRbSMmuj/csFOYazfqu57hJ32CBzdpN3phuOkofCvKV+SSBLh0BwShm0Yf9kQ2WzsTzFPKU9jagA6oT2QdXAgJ6GILcS0mujUkur0BeYXa76FqqmUat9FosLq6SqVSoV6vt2wn2846uA0stXuqzbVrG1rie0Li2QK12wtUThSotbcbLbW1G53debvRjQzbNPpOqBLENdsS1G1s8UrD9+Fjr7JVGmUdXAgJ6GIT6XpsuBhCBEQW7osIzrOsx+vUbqoRLLi9++mb6m62k8WLETnlktdq9/vU7i5S/VqBcLF1Gtt/WNgs7lK8KEzaje5dmsyWTs3Oz88394OPGu/iPHap4Uq/euNVAlYamwixudH/Xy56Jl2PDaOIiq1QYZ2GbRCeisnP5vGsR87b27RmuKJYv6lE4+4ylbvz2EZmr68fUz56Zio9P9+9dqNpNrrv+0xNTTWz0Uc9OHglb2xKwGaz0dMALtPoQmxsogN6XI0JTzZQdU+292SkPcJXHjxNvVEniAOUVeRUDpRH3npuW9cuBrBbthvdF1K+qE7xwjqlSwPyF3TnTzQdhafJbPv372d2dnbXSXnDbFRLwEpddCH2ZvzezbYprsaEn60TRzG53ORu74lrMfU7qqyurFNVFcL5iBBXDCWd0vanfexp6xpuA3hsupWr3abtRn1L6dKAqYc18C+o4s/FZ7aKnbu3P88wDLHWUiwWKZVKzM/PS1GXIZNtdiR10YXYm4kN6NHtwcRu7wmCgNXVVdaW16jeUSWyITnl3kTVKuQP51F5hV1wyWnWAzWvsBULEagFhXcwv2GxlWa70RMFKic3aDd6havOlrYbnfLLrN4Td9wqtl3ZtfBSqTTWo/BR1V5WVbaTCdE9E/tOZ6tJMM8sy47j9h5rLUEQtHQlS5PYOGVRCnx1Zu7cKuuS1M7Ln70ne3rjQNvSbvRkgXh9g3ajxxv4B85eC1cFr7lVbCfSfeHpWniakS4BYjjIdjIh+mdiA7oqQ1xp3eo0Dtt7svvAa7UajUZjwyz0KAhbSpzC2XXLs3uyWx8naTd6okD1hE/97v61G82WV52bm2N+fn4kM9LHlUyjCzEYExvQc5f5xKcibGwhx8hu74njmFqt1iyl2l7IZdNtZDusW97SbvREgWi1td1o6dKgmZHunxvRrffvbEJbqVTi4MGDI7MvfBLINLoQw2G0olcXeWWP/KOKeDfXoK5GZntPFEVUKhXW1taaLUWzpVR30tBkq7rl1kJwf65ZI/2sdqNzcbNTWflIgFfqTXGXqakp6VI2RDabRp+amqJSqQz6EIWYSBMb0CEJ6scKQ72eFwQBp0+fZm1tjUaj0Wxmkk5h7iXhq1PdcqZzVG9x1dmqJ3yi05kPB96ZdqNTxwP8h3VvFA4ylT7MZBpdiOE30QF92KQJbGkZ1UajQbFYpFKpdCWAd5RXxF6Byp0ugNfu6NBu9JgL4B3bje5BNiu9UChw3nnnMTMzI1PpQ6C9qItMowsx/CSgD1C6/r22ttYxgS1tJ9rtIN7SbvREgXCxtd1o8cKkRvoVDQrnb95udKey6+Hlcpn5+XmmpqY4cOBAX3qFi87ae4RLURchRo8E9D7Krn83Go1mK9FtJbDtUbDoUb3ZFXap3eZjw8wofKqt3eh090bhcGZrWZ48xdUC56qDlGaK5M7xt8xZ2G6XN7Fz0iNciPEiAb1H0pHo2toa6+vrBEFAEARd6QW+rcfPtButnCwQPtg6jV24INNu9PDe2o12kia1pS1H58vzeF+2zWI+dtkSbVGZL20Ok/7OpFbz6xapjS7EeJOA3iXt0+dBEDTXh9M14V5XLAuXPaonXADvZ7vR5uNnktrm5+eZn59vnnN4okGs7I4q86XNYSaxml+3tLcY9X1fchSEGFMS0HcpDEOq1epZ28fS6fN0/buXbAT1u/LNjPTg/taXs1ftRrPC0NV9TxuezM3NdQwYtkZLj27YujJf++/EYYxdjokXkyAl0+9nSbeUpcs3aY9wmUYXYvxJQN+GbPZ5tVpt2T7Wj+nzrHBFUT3pAnj1Fh9bzzQ6KVjKR5N94V1uN5pqr5d+8OBBpqent5y2VSWIa7YlQG9VmS/7O3EYY++LiLF4RY/gtgbqRIB3WY78kcmdgpfSqkKIlAT0DqIoolqtsr6+Tq1Wo16pEz0U4kU5vIKrZ96v/dE2gtqdmXaj97WNwg+GlK9ISqxeHKJ68Iq2V2o799xztxXEs7yL80QP1olWQlSssJ7Fm8ttWpnPuziPXWpgY+tG5liUVdhaDJ4rhBPfHRGtTda6ensym+wJF0KABHTAFW/ZKHmNENS91nUjQ2GrrgOZd3jjbmN7Fa0pqiddAK/dWiCqlJrXKd9SujxTYnX/9kbhNrDNAjL4bNnNrH172b59+7rQ9ORMmVn37+br+F7Jg6sLxHeFxIsxXtlzpXob6sxxxK6E/Divq7ePwiWZTQjRyUQHdGstt9xyC0DLtHl2+jxaTBKzMoEo25GsK8cRQ+PreSrpKPzrudZ2o+dEzepsxUtdu9Hm724jUNsgaYOanIcNO38oSYN4Pp9vBvFyudyVkV98V4jyPbyDrVPuWwVir+Q1r4+XY+IHI9LS89ZaVE6NZZe89sps6Vq4EEJsZOIDOrD59HnAlh3JdiOqJKPwk25veEu70bylfJlbC1+4Ok84vdb5+LcZqOPFaMMPJd65uWYQ72XN9N0kxWU1p9892xzYKwtqnzc2XfKkMpsQYi8mOqBvyw47km1kW+1G07Xwy860Gy1OTxOud77PzQJ1y+xB24cSiyWyEblGjqmpKfbv30+pVGq/+67aTVJcVjr9zq0e0e0NyCvUfg/lqZHskgetW8qkMpsQYq9G712wz7bqSLaZuKao3uI318Nb2o3mLKVL3Dp4+XiAf3AXjU62O3vgQxxYYiJyKkc5X2Zffp7yOWXy5xeIazHhiUZPq7FlE9zS/eQ7DcReycO7qkDuSH4kq8e1r4XLljIhRDdJQN9Cp45kGyWUZduNVk4UqN/V3m40SgJ4w7UbLe7x4DaZPbCBJVoMiWoxOT9HmSL7y/sp+aWWYNqvamzZBDdbY0/tarPr6sMuOwpPW8HKKFwI0QsTG9Djakxwa53wvhCvuHGQBhfEN0qAi+tQu80fSLvRTrMHxBBNxdh7LOVciX2l/RS9ggv0cwqi1t7v4YlG36qxjVIg3q1sYZe0vGo6Ci+VSqytdc6HEEKIvZrIgB5XY8LP1omthTDJtt7mVjRrITzlNQN4p3ajU8mWstKR7rYbbZfOHkSLIWEtJFfKMXNoltnVGYozhdYkNB883yN/VWtA3Wuy2qRr71ImhV2EEIMykQE9uj1wo1K1va1ocSNtN+qCeLjU1m70oqA5ld7tdqMbsdYSRRG5fI7ZS+bYv38/xWIRpRTBFxrgtX6Q2ChI7zVZbRK1j8KlS5kQYhgMNKBrrXPA24CrgBzwM8aYL/b6cW01GZXGZ4JeezJZ8JDXDOC12/vbbnTD484Ue0m3mHXaJ76TIN2NZLVxJ+VVhRCjYNDv2s8EPGPM92qt/zvwFuCJvX5QVYa42hqE49DSeKBI/fMlKicKhKfa2o0eDikfc2vhhR60G91Itnb6diu27SRIdzNZbZy0j8KlvKoQYtipdP1vELTWbwU+TjI6B640xlzcdpvrgOsAjDGPaDQazevy+XyzZedORJWQlU8us7oc8J+fv43G3WWq97SNwsuWmW+Kmb0yYuaKmPzcbs5wd7Ij8ZmZGRYWFpiZmdlRMImqIbVbq8QVizelKB0pkysP+vPbxnb7WnZLGsDTY+nFvvBBn2M/TMI5wmSc5yScI4zGeRYKBYAtA8AwvMO/Cvgo8FTg5vYrjTHvAt6V/GhPnTrVvO7AgQNkf96Ju0+tc9PbV4D55mX++WEzoa144Zl2o3WgvkFxl25KW5GWy2Xm5+eZmZnhnHPOYXFxkaWlpZ3f4QVnvq1XA6h271i7bWFhgcXFxb4+5kZr4WknvW7by9/rqJiEc4TJOM9JOEcYjfM8dOjQtm436ID+OcA3xrxGa/2dwFf79cDnfEuRXFnhX1pl+oqQ8rHetBvdShRFADtqRSp2r31fuKyFCyHGxaAD+o3AE7XWnwRC4EX9euDpQ3me8KcH+NJXvti3VqipdHqnWCxyzjnnMDs7K0GlRzbbFy6EEONkoAHdGBMA1w7q8b18/97UoyjCWkuxWGT//v3Mzc1J96weSfuF53I5GYULISbGoEfoYy2OY6IoolQqMT8/z/z8fEtrVtEdaXEXqZEuhJhkEl26LA3ihUKBubk5FhYWJIj3gHQqE0KIVhJpuiDdZub7PnNzrmpbv9flx530CxdCiM1JQN+lNIjn8/m+9RSfNOkoXCnVTGiTUbgQQnQmAX0HsgVfyuUy+/fv71h6VexOtsRqGsB935fnVwghtkEC+jZkC76ce+65O67aJjaWBvE0Kz2dShdCCLEz8s65gXSveKlUkr3iXdapuMvs7Cz1en3ARyaEEKNLAnqG7BXvjey2MinuIoQQvTHxAT2OY4IgoFgsMjc3x759+2TKtwtkW5kQQvTXREcupRQLCwssLCzINrM9ak9ok21lQgjRXxMf0A8ePCgjx13q1DNcZjeEEGIw5N1X7Ih0KxNCiOEkAV1sSuqkCyHEaJCALs4iFdqEEGL0SEAXLQltnudJhTYhhBhBEtAnVPve8GKxKFPpQggxwiSgTxDZGy6EEONLAvoYk2YnQggxOSSgj5lOe8NlKl0IIcafBPQxIFPpQgghJKCPIJlKF0II0U4C+oiQrHQhhBCbkYA+xOI4JooiAJlKF0IIsSkJ6EMmiqKWAi9zc3MEQTDowxJCCDHkJKAPWPtUeqFQaJlKlyl1IYQQ2yEBfQDat5ZJrXQhhBB7JQG9T6TtqBBCiF6SgN4j6VQ6QC6Xk7ajQgghekoCehe1T6Xn83nyeXmKhRBC9J5Emz2SqXQhhBDDQAL6LmSDuFRpE0IIMQwkoG/DVlvLhBBCiEGTgL4B2VomhBBilEhAz5D1cCGEEKNq4gN6GsSla5kQQohRNtEBXSnF9PQ0nudJEBdCCDHSJjqggxuZCyGEEKNuoAFda30DcDWwDJwwxvz0II9HCCGEGFUqLU86CElAv8EY84+b3OY64DoAY8wjGo1G87p8Pk8Yhj0+ysGahHOEyThPOcfxMQnnOQnnCKNxnoVCAWDLdeG+BHSt9TuBR7Rd/FPAc4DvBqrAO40xN25xV/bee+9t/nDgwAFOnTrVzUMdOpNwjjAZ5ynnOD4m4Twn4RxhNM7z0KFDsI2A3pcp902m0r8AoLXeB/yn1vr/GGPW+3FMQgghxDgZlk3WEVAHGlvdUAghhBBnG3RS3FuBh+MC+kuMMcEgj0cIIYQYVQMN6MaYlw3y8YUQQohxMdAs910YqYMVQgghumTLpLhhWUPfLpX90lp/rv2ycfuahHOclPOUcxyfr0k4z0k4xxE7zy2NWkAXQgghRAcS0IUQQogxMOoB/V2DPoA+mIRzhMk4TznH8TEJ5zkJ5whjdJ6jlhQnhBBCiA5GfYQuhBBCCCSgCyGEEGNBAroQQggxBgZaKW43tNY54G3AVUAO+BljzBcHe1TdpbWuAZ9Jfnwr8HXgzbjz/VdjzC8N6th2S2tdAl4I/CJwvTHmBq31t9N2XqP8+m5wjtcDzwK+AZw2xvzwKJ8jgNZ6CpdIdAEwBTwXmGe8XstO5/hsxuy1BNBafwwoACVcq+oyY/RaQsdzfAZj+FqO4gj9mYBnjPle4NeAtwz4eHrhG8aY70u+PoR7Y7nWGPNdwNVa68cM+Ph24zxcm9w/y1zW6bxG+fXtdI4Ab0heyx9Ofh7lc8QYUwF+0xjzWODdwIsZs9dyg3OEMXstAYwx1yTn+VHgMYzZawkdzxHG8LUcuSz3pKHLx0k+QQFXGmMuHuxRdZfW+lbgXuAU8NvAO4BrgDcB3wzcaIx56+COcPeSEeudwF/hXseW8wIuY8Rf3/QckxH6y4FrcYH+/caYt43T33ByrhZ4CmP4WkLzHJdw1brG7rXUWj8ZeDWwD/he4MOM2WvZdo7fA2jG8LUcxRE6wKuAq4GnAuGAj6XrjDFHjDHfA7wPeC1wKPn+bZw9+htlG53X2Ly+xpjfN8Y8Evh+4Fqt9fHkqpE/R63104BvBd7OmL6WmXP8w3F9LY0xH0nO6+dwg4exey3bz3FcX8tRDOifA/7LGPMa4JHAVwd8PL2kgDuANeAFwE3Ak4BPD/KgusEYs0Tn8xrX1zdK/l1lDM5Ra/1TwI8CP2qMeZAxfC3bzjH75j5Wr2XGMlBhDF/LjGXcqDw1Vq/lKE65+8B7gAtxn6BeZIy5ZbBH1T1a63OBDwB14AHgJbgRwuuBAPh48gc3UrTWh3FTeYdw5/bPwA20ndcov74bnONXgSfipqXfYYwxo3yOAFrrRwL/BnwKiIEG8AbG67XsdI7/wPi9lpfi/h/WgBXg5cBxxuu17HSOz2XMXksYwYAuhBBCiLON4pS7EEIIIdpIQBdCCCHGgAR0IYQQYgxIQBdCCCHGgAR0IYQQYgxIQBdiDGmtD2qtveT7b0v+PaS1VoM9MiFEr0hAF2LMJLUMPoEr2wnwfq11AXgO8BGtdTlz26u01h/O/Dyvtf4MW9Ba+1rrf9ZafznZv7vR7Z6mtX7Crk9GCLFtsg9diDGitT6Aa0Dxe8aYP08u+wvgr4GTuD+XL1AAAARGSURBVFKXXwOux9UnL+E6it2W3MWHge8Afj5ztw8llf3Sx/CAP8bVxV7GDQx+yhgTZ27zIlyt7FlcoY6vAA8BPq7wziW4wkkV4DPGmBd16SkQYmLJCF2IMZG0o/0UcBHwoNb6P7XWn8MF7HcDz8e1xjwfeH1Sy/pZwCeNMY9Mfp4GZnA9BF4LfBDXEjZ9jDLwp7hqYtcCPw2cC/x50nI0dT6uxvufGGMWjDH/HTiMa4rxBOAuXFnRJwAv68HTIcTEkRG6EGNCa30drlzpW4DnGWPuyVz3Dlw52l/BBeP/Am7Bfagv4MpifgJ4HPAo4C+NMddorf8O+HljzBe11lfhmnXcDvwPY0w1uW8fN2L/TuA6Y8wnM4/5FFx3PXCNL16AK2Wc1tM+D3iFMeZ9XX9ChJgw+UEfgBCiO4wx7wLQWn8BF7yzXoybZk+vexlwc9ttGsA9uOYcR7TWOeBzSTD/FeCluKn7JwFf1VrP4D4MLCa/fyPwLq31B3EzBZcAP5jc31/hmn0A/DjwMNwHg6fu9byFEI4EdCHGz1XAP2itO11XAH4D+Ajw9LbrThhjXqe13oebti8YY341ue5/A39gjKmkN9ZavwC4yhjz8sxlr8JN2X8PrgnGu5Ljea0xZjFzTE8A/n5PZymEaCEBXYjxc9wYcwmA1voi4Kgx5u+Tn2/AJaO9EZf81kJr/TxgHtdp61HAPyVXPQZ4U9uHhCnA11o/s+1ufhm4G/g14EHcKPyFWusPAevJbYq4aX4hRJdIUpwQ4+0q4GntFxpjXgL8L+AngVO4UfvjjTFfwk2Tvw43NZ7e/kZjzGHgcuC7k+9/BXh38v1jjTGHk68bcSP8b8Il5L0Zt27/q5xxLi7rXQjRJRLQhRhvjwJKWuti9kKt9dW4LPUguegg8Fta63OAa4DfAa5ORvhZrwPe1HZf88A/aa1bPgAAgTHm0caYR+PW6z3cervFjfhDYJUzo3YhxB5IQBdi/DxLa32O1vqVwLNxe72/pLV+J/CXuEz4F+HWuGeBCLcV7R+B9+K2tAW4kfUNWusSgNb6KcCP4RLsmowxp3GJcm/TWj9/i2N7GfA84ATwfuB+Y8yH9nrCQggJ6EKMjaR627tx+78/jVunfrQx5mW4qfdP4Na3PwX8OvAaXHLcR5KiMD8EfNkY80EAY8z7cQVn3qC1/g7gBuCZxpgHkoc8iMtgxxjzX7ig/hat9ROT63Na688kleeejKtc9zFc0tyTgf+O+/DxWz16SoSYKJIUJ8SYMMYEWus/Ah4wxtzddl0dl6n+v7XWlxpjHgSe0XYXrzTG3Nd22U/jkt9KgDbGfEZr/Wbc+noN9yEgfYx/01o/HvhqctGTjTHplja01hcCtxljPp9cVMcF9Et2f9ZCiJQUlhFCCCHGgEy5CyGEEGNAAroQQggxBiSgCyGEEGNAAroQQggxBiSgCyGEEGPg/wcou/j/Hb7jKAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dalta_pred = t_half_alpha * s_err * np.sqrt(1 + 1.0/n + np.power(xs_pred - x_bar, 2)/x_var)\n",
    "\n",
    "lower_95_percent = ys_pred - dalta_pred \n",
    "upper_95_percent = ys_pred + dalta_pred\n",
    "\n",
    "ax.fill_between(xs_pred, lower_95_percent, upper_95_percent, color='#888888', alpha=0.1)\n",
    "fig"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
